I have a program that works in VS C++ and does not work with g++. Here is the code:
#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>
#include <vector>
#include <cstdio>
#include <algorithm>
#include <set>
#define EP 1e-10
using namespace std;
typedef pair<long long, long long> ii;
typedef pair<bool, int> bi;
typedef vector<ii> vii;
// Returns the orientation of three points in 2D space
int orient2D(ii pt0, ii pt1, ii pt2)
{
long long result = (pt1.first - pt0.first)*(pt2.second - pt0.second)
- (pt1.second - pt0.second)*(pt2.first - pt0.first);
return result == 0 ? 0 : result < 0 ? -1 : 1;
}
// Returns the angle derived from law of cosines center-pt1-pt2.
// Defined to be negative if pt2 is to the right of segment pt1 to center
double angle(ii center, ii pt1, ii pt2)
{
double aS = pow(center.first - pt1.first, 2) + pow(center.second - pt1.second, 2);
double bS = pow(pt2.first - pt1.first, 2) + pow(pt2.second - pt1.second, 2);
double cS = pow(center.first - pt2.first, 2) + pow(center.second - pt2.second, 2);
/* long long aS = (center.first - pt1.first)*(center.first - pt1.first) + (center.second - pt1.second)*(center.second - pt1.second);
long long bS = (pt2.first - pt1.first)*(pt2.first - pt1.first) + (pt2.second - pt1.second)*(pt2.second - pt1.second);
long long cS = (center.first - pt2.first)*(center.first - pt2.first) + (center.second - pt2.second)*(center.second - pt2.second);*/
int sign = orient2D(pt1, center, pt2);
return sign == 0 ? 0 : sign * acos((aS + bS - cS) / ((sqrt(aS) * sqrt(bS) * 2)));
}
// Computes the average point of the set of points
ii centroid(vii &pts)
{
ii center(0, 0);
for (int i = 0; i < pts.size(); ++i)
{
center.first += pts[i].first;
center.second += pts[i].second;
}
center.first /= pts.size();
center.second /= pts.size();
return center;
}
// Uses monotone chain to convert a set of points into a convex hull, ordered counter-clockwise
vii convexHull(vii &pts)
{
sort(pts.begin(), pts.end());
vii up, dn;
for (int i = 0; i < pts.size(); ++i)
{
while (up.size() > 1 && orient2D(up[up.size()-2], up[up.size()-1], pts[i]) >= 0)
up.pop_back();
while (dn.size() > 1 && orient2D(dn[dn.size()-2], dn[dn.size()-1], pts[i]) <= 0)
dn.pop_back();
up.push_back(pts[i]);
dn.push_back(pts[i]);
}
for (int i = up.size()-2; i > 0; --i)
{
dn.push_back(up[i]);
}
return dn;
}
// Tests if a point is critical on the polygon, i.e. if angle center-qpt-polygon[i]
// is larger (smaller) than center-qpt-polygon[i-1] and center-qpt-polygon[i+1].
// This is true iff qpt-polygon[i]-polygon[i+1] and qpt-polygon[i]-polygon[i-1]
// are both left turns (min) or right turns (max)
bool isCritical(vii &polygon, bool mx, int i, ii qpt, ii center)
{
int ip1 = (i + 1) % polygon.size();
int im1 = (i + polygon.size() - 1) % polygon.size();
int p1sign = orient2D(qpt, polygon[i], polygon[ip1]);
int m1sign = orient2D(qpt, polygon[i], polygon[im1]);
if (p1sign == 0 && m1sign == 0)
{
return false;
}
if (mx)
{
return p1sign <= 0 && m1sign <= 0;
}
else
{
return p1sign >= 0 && m1sign >= 0;
}
}
// Conducts modified binary search on the polygon to find tangent lines in O(log n) time.
// This is equivalent to finding a max or min in a "parabola" that is rotated and discrete.
// Vanilla binary search does not work and neither does vanilla ternary search. However, using
// the fact that there is only a single max and min, we can use the slopes of the points at start
// and mid, as well as their values when compared to each other, to determine if the max or min is
// in the left or right section
bi find_tangent(vii &polygon, bool mx, ii qpt, int start, int end, ii center)
{
// When query is small enough, iterate the points. This avoids more complicated code dealing with the cases not possible as
// long as left and right are at least one point apart. This does not affect the asymptotic runtime.
if (end - start <= 4)
{
for (int i = start; i < end; ++i)
{
if (isCritical(polygon, mx, i, qpt, center))
{
return bi(true, i);
}
}
return bi(false, -1);
}
int mid = (start + end) / 2;
// use modulo to wrap around the polygon
int startm1 = (start + polygon.size() - 1) % polygon.size();
int midm1 = (mid + polygon.size() - 1) % polygon.size();
// left and right angles
double startA = angle(center, qpt, polygon[start]);
double midA = angle(center, qpt, polygon[mid]);
// minus 1 angles, to determine slope
double startm1A = angle(center, qpt, polygon[startm1]);
double midm1A = angle(center, qpt, polygon[midm1]);
int startSign = abs(startm1A - startA) < EP ? 0 : (startm1A < startA ? 1 : -1);
int midSign = abs(midm1A - midA) < EP ? 0 : (midm1A < midA ? 1 : -1);
bool left = true;
// naively 27 cases: left and left angles can be <, ==, or >,
// slopes can be -, 0, or +, and each left and left has slopes,
// 3 * 3 * 3 = 27. Some cases are impossible, so here are the remaining 18.
if (abs(startA - midA) < EP)
{
if (startSign == -1)
{
left = !mx;
}
else
{
left = mx;
}
}
else if (startA < midA)
{
if (startSign == 1)
{
if (midSign == 1)
{
left = false;
}
else if (midSign == -1)
{
left = mx;
}
else
{
left = false;
}
}
else if (startSign == -1)
{
if (midSign == -1)
{
left = true;
}
else if (midSign == 1)
{
left = !mx;
}
else
{
left = true;
}
}
else
{
if (midSign == -1)
{
left = false;
}
else
{
left = true;
}
}
}
else
{
if (startSign == 1)
{
if (midSign == 1)
{
left = true;
}
else if (midSign == -1)
{
left = mx;
}
else
{
left = true;
}
}
else if (startSign == -1)
{
if (midSign == -1)
{
left = false;
}
else if (midSign == 1)
{
left = !mx;
}
else
{
left = false;
}
}
else
{
if (midSign == 1)
{
left = true;
}
else
{
left = false;
}
}
}
if (left)
{
return find_tangent(polygon, mx, qpt, start, mid+1, center);
}
else
{
return find_tangent(polygon, mx, qpt, mid, end, center);
}
}
int main(){
int n, m;
cin >> n >> m;
vii rawPoints(n);
for (int i = 0; i < n; ++i)
{
cin >> rawPoints[i].first >> rawPoints[i].second;
}
vii polygon = convexHull(rawPoints);
set<ii> points(polygon.begin(), polygon.end());
ii center = centroid(polygon);
for (int i = 0; i < m; ++i)
{
ii pt;
cin >> pt.first >> pt.second;
bi top = find_tangent(polygon, true, pt, 0, polygon.size(), center);
bi bot = find_tangent(polygon, false, pt, 0, polygon.size(), center);
// a query point is inside if it is collinear with its max (top) and min (bot) angled points, it is a polygon point, or if none of the points are critical
if (!top.first || orient2D(polygon[top.second], pt, polygon[bot.second]) == 0 || points.count(pt))
{
cout << "INSIDE" << endl;
}
else
{
cout << polygon[top.second].first << " " << polygon[top.second].second << " " << polygon[bot.second].first << " " << polygon[bot.second].second << endl;
}
}
}
My suspicion is there's something wrong with the angle function. I have narrowed it down to either that or find_tangent. I also see different results in g++ when I switch from double to long long in the angle function. The double results are closer to correct, but I can't see why it should be any different. The values I'm feeding in are small and no overflow/ rounding should be causing issues. I have also seen differences in doing pow(x, 2) or x*x when I assign to a double. I don't understand why this would make a difference.
Any help would be appreciated!
EDIT: Here is the input file: https://github.com/brycesandlund/Coursework/blob/master/Java/PrintPoints/points.txt
Here is the correct result: https://github.com/brycesandlund/Coursework/blob/master/CompGeo/CompGeo/correct.txt
Here is the incorrect result: https://github.com/brycesandlund/Coursework/blob/master/CompGeo/CompGeo/fast.txt
