iv. Polygons

Polygon Area

For a simple polygon with vertices $p_0,\dots ,p_{n-1}$ in order:

$2A={\sum }_{i=0}^{n-1}(x_i{y}_{i+1}-x_{i+1}{y}_i)={\sum }_i\mathrm{cross}(p_i,p_{i+1}),$

with indices mod $n$. The signed value is positive for counter-clockwise order.

ll area2(const vector<P>& p) {                        // TWICE the signed area
    int n = p.size();
    ll s = 0;
    for (int i = 0; i < n; i++) s += cross(p[i], p[(i + 1) % n]);
    return s;                                          // > 0 means CCW
}
double area(const vector<P>& p) { return fabs((double)area2(p)) / 2.0; }
 
PD centroid(const vector<P>& p) {                      // area-weighted centroid of the region
    int n = p.size();
    double cx = 0, cy = 0, a2 = area2(p);
    for (int i = 0; i < n; i++) {
        ll c = cross(p[i], p[(i + 1) % n]);
        cx += (double)(p[i].x + p[(i + 1) % n].x) * c;
        cy += (double)(p[i].y + p[(i + 1) % n].y) * c;
    }
    return {cx / (3 * a2), cy / (3 * a2)};
}

Keep twice the area

$2A$ is always an integer for integer coordinates. Work with it throughout, halve only when printing, and you never touch floating point.

Proof: the shoelace formula and the area-weighted centroid

Area. Pick any point $O$ (the origin, say) and sum the signed area of triangle $O{p}_i{p}_{i+1}$ over all edges: ${\sum }_i\frac{1}{2}\mathrm{cross}(p_i-O,p_{i+1}-O)$. For an edge of the polygon that is “outer” relative to $O$ (the triangle $O{p}_i{p}_{i+1}$ lies outside the polygon, contributing area that should be subtracted), traversing the polygon CCW makes that triangle’s own orientation CW, so its signed area is already negative — it cancels automatically against the true polygon area, which is why the formula is origin-independent and correct for any simple polygon, convex or not, and for any choice of $O$ including $O=$ origin (giving the stated formula with $p_i,p_{i+1}$ directly, since $\mathrm{cross}(p_i-O,p_{i+1}-O)=\mathrm{cross}(p_i,p_{i+1})$ when $O=0$). A full proof by induction on $n$: triangulate off one ear $p_{k-1}{p}_k{p}_{k+1}$ (a valid ear exists by the two-ears theorem, see below); the shoelace sum over the $(n-1)$-gon with $p_k$ removed plus the signed area of the ear telescopes (the two new edges $p_{k-1}{p}_{k+1}$ appear with opposite signs in the two sums and cancel) to the $n$-gon’s shoelace sum, while the true area also decomposes as (area of smaller polygon) + (area of the ear) since the ear ray is cut off cleanly — so if the formula holds for $n-1$ vertices it holds for $n$; the base case $n=3$ is the signed-triangle-area identity from Primitives.

Centroid. The centroid of a region is the area-weighted average of the centroids of any triangulation: $\bar{p}=\frac{{\sum }_T\mathrm{area}(T){\bar{p}}_T}{{\sum }_T\mathrm{area}(T)}$. Triangulate as the fan $\{O,p_i,p_{i+1}{\}}_i$ from the origin (signed, so pieces outside the polygon again cancel exactly as in the area proof); triangle $O{p}_i{p}_{i+1}$ has centroid $\frac{1}{3}(p_i+p_{i+1})$ (origin contributes $0$) and signed area $\frac{1}{2}\mathrm{cross}(p_i,p_{i+1})$. Summing signed-area-weighted centroids and dividing by the total signed area $A$ gives $\bar{p}=\frac{{\sum }_i\frac{1}{3}(p_i+p_{i+1})\cdot \frac{1}{2}\mathrm{cross}(p_i,p_{i+1})}{A}=\frac{{\sum }_i(p_i+p_{i+1})\mathrm{cross}(p_i,p_{i+1})}{6A}$, the formula above. $■$

What the sign gives you

• $2A>0$ $\to$ counter-clockwise; $2A<0$ $\to$ clockwise; $2A=0$ $\to$ degenerate.
• Normalise once at input: if (area2(p) < 0) reverse(p.begin(), p.end());

Many algorithms (point-in-polygon, hull merging, Minkowski sums) silently assume CCW.

Requirements

The polygon must be simple — no self-intersections. For a self-intersecting polygon the formula returns a signed sum with regions counted by winding number, which is occasionally what you want but usually is not.

Related quantities

Quantity	Formula
Perimeter	$\sum |p_{i+1}-p_i|$
Centroid	$\frac{1}{6A}\sum (p_i+p_{i+1})\mathrm{cross}(p_i,p_{i+1})$
Pick’s theorem	$A=I+\frac{B}{2}-1$ for a lattice polygon
Boundary lattice points	$\sum \gcd (|\Delta x_i|,|\Delta y_i|)$
Interior lattice points	$I=A-\frac{B}{2}+1$
Is it convex?	all orientations agree
Bounding box	min/max of the coordinates

The centroid formula is the area-weighted one (the centre of mass of the region), not the average of the vertices. The two differ for non-uniform vertex spacing, and problems usually mean the former.

Triangle area

$2A=|\mathrm{cross}(b-a,c-a)|.$

Also useful: Heron’s formula $A=\sqrt{s(s-a)(s-b)(s-c)}$ with $s$ the semiperimeter — but it is numerically unstable for thin triangles. Prefer the cross product when coordinates are known.

Area of a union of shapes

Shapes	Method
Union of rectangles	sweep line + segment tree with counts, $O(n\log n)$
Union of circles	angular sweep on each circle’s boundary, $O(n^2\log n)$ — 13-circle-union-and-angular-sweep
Union of general polygons	polygon clipping (Weiler-Atherton, or a Boost/Clipper library) — 11-polygon-clipping-and-union
Intersection of two convex polygons	$O(n+m)$ by simultaneous traversal, or half-plane intersection — 09-half-plane-intersection
Intersection of a convex polygon and a half-plane	Sutherland-Hodgman clipping, $O(n)$
Circle-polygon intersection area	sum the signed areas of circle-triangle pieces per edge

Sutherland-Hodgman clipping

vector<PD> clipLeft(const vector<PD>& poly, P a, P b) {  // keep the side with orient(a,b,*) >= 0
    vector<PD> res;
    int n = poly.size();
    for (int i = 0; i < n; i++) {
        PD cur = poly[i], nxt = poly[(i + 1) % n];
        double s1 = cross(toPD(b - a), cur - toPD(a)), s2 = cross(toPD(b - a), nxt - toPD(a));
        bool in1 = s1 >= 0, in2 = s2 >= 0;
        if (in1) res.push_back(cur);
        if (in1 != in2) {                              // interpolate the crossing point
            double t = s1 / (s1 - s2);
            res.push_back(cur + (nxt - cur) * t);
        }
    }
    return res;
}

Clipping a convex polygon by each of another’s edges gives their intersection in $O(nm)$ — enough for most problems; the full clipping/union toolkit (including non-convex clip windows) is its own note.

Pick’s Theorem

For a simple polygon whose vertices are lattice points:

$A=I+\frac{B}{2}-1,$

where $I$ is the number of interior lattice points and $B$ the number of lattice points on the boundary.

Computing the three quantities

Area — shoelace, exactly (see above). Boundary points — each edge from $a$ to $b$ contributes $\gcd (|\Delta x|,|\Delta y|)$ lattice points, counting one endpoint per edge:

ll boundaryPoints(const vector<P>& p) {
    int n = p.size(); ll B = 0;
    for (int i = 0; i < n; i++) { P d = p[(i + 1) % n] - p[i]; B += __gcd(llabs(d.x), llabs(d.y)); }
    return B;
}

Interior points — rearrange Pick: $I=A-\frac{B}{2}+1=\frac{2A-B+2}{2}$, computed from area2 and B with no floating point:

ll interiorPoints(ll twiceArea, ll B) { return (llabs(twiceArea) - B + 2) / 2; }

Proof: the boundary gcd count

The lattice points strictly between $(0,0)$ and $(\Delta x,\Delta y)$ are exactly the integer multiples $k\cdot (\frac{\Delta x}{g},\frac{\Delta y}{g})$ for $k=1,\dots ,g-1$, where $g=\gcd (|\Delta x|,|\Delta y|)$ and $(\Delta x/g,\Delta y/g)$ is the primitive direction (any lattice point on the segment has coordinates that are both multiples of $\Delta x/g$-scale steps simultaneously in both axes exactly when the step divides both $\Delta x$ and $\Delta y$, and the largest such step is $1/g$ of the vector by definition of gcd). That is $g-1$ interior points, plus counting one endpoint ($a$, say) per edge gives $g$ points charged to that edge; summing over all edges charges every boundary lattice point exactly once (each vertex is the “charged” endpoint of exactly one edge going around). $■$

Proof: Pick's theorem

Define the defect $D(P)=A(P)-I(P)-\frac{B(P)}{2}+1$ for a simple lattice polygon $P$; the claim is $D(P)=0$ always.

Base case. A primitive triangle (lattice vertices, no other lattice point on its boundary or interior) has area exactly $\frac{1}{2}$: by a unimodular integer change of coordinates (translate a vertex to the origin; the other two edge vectors $u,v$ from it satisfy $\gcd$-primitivity on each edge and $|\mathrm{cross}(u,v)|=1$ — a classical fact: if $|\mathrm{cross}(u,v)|>1$ the triangle contains an extra lattice point, found via the extended Euclidean algorithm on $u,v$ as a basis change), the area $\frac{1}{2}|\mathrm{cross}(u,v)|=\frac{1}{2}$. Then $A=\frac{1}{2}$, $I=0$, $B=3$, so $D=\frac{1}{2}-0-\frac{3}{2}+1=0$.

Additivity. If a simple lattice polygon $P$ is split by a diagonal into lattice polygons $P_1,P_2$ sharing an edge with $k\ge 2$ lattice points on it ($k-2$ of them strictly between the endpoints): $A(P)=A(P_1)+A(P_2)$; every lattice point of $P$ is interior/boundary to $P_1$ or $P_2$, and the $k-2$ points strictly inside the shared edge move from “boundary of $P$” to “interior of both pieces, counted twice”, while the two shared endpoints stay boundary points of $P$ but become boundary points of both pieces. Careful bookkeeping ($I(P_1)+I(P_2)=I(P)+(k-2)$, since those interior-to-the-cut points were boundary of $P$ and are now interior to the union of pieces counted once each in $I(P_1)+I(P_2)$… after the standard accounting) gives $D(P)=D(P_1)+D(P_2)-1$ is the naive guess, but the correct invariant to track is $A-I-B/2$ (dropping the $+1$), which is additive: $B(P_1)+B(P_2)=B(P)+2k-2$ (the shared edge’s $k$ points are counted in both pieces’ boundary, i.e. $2k$ times, vs. once each as boundary-or-now-interior in $P$ contributing differently)… the clean way to see it without heavy casework: both $(A-I-B/2)$ for $P$ and the sum for $P_1,P_2$ change consistently under adding one triangle at a time (ear clipping), because adding a single ear (triangle with two old vertices and the new structure sharing exactly one full edge with the rest, no new vertices introduced beyond the original polygon’s) changes $A$ by the ear’s area, changes $B$ only by points on the newly exposed diagonal that were already lattice points of the original polygon (none are new, since we never introduce Steiner points — vertices stay the original lattice points) and swaps some previously-boundary points of the whole polygon to interior/shared status whose net effect on $I+B/2$ matches the ear triangle’s own $I+B/2$ contribution exactly when the ear is primitive. Concretely: triangulate $P$ using only its own vertices as triangle vertices (always possible for a simple polygon — the two-ears theorem) into $n-2$ triangles, then subdivide any non-primitive triangle further using its own interior/boundary lattice points as additional vertices (always possible: a non-primitive lattice triangle has an extra lattice point on its boundary or inside, and connecting it splits the triangle into strictly smaller-area lattice triangles; recurse, terminating since area strictly decreases and is bounded below by $\frac{1}{2}$ for a non-degenerate lattice triangle). Now $P$ is triangulated into $N$ primitive triangles for some $N$. Each primitive triangle has $D=0$. Gluing primitive triangles back together one shared edge at a time, each gluing step merges two pieces along a lattice segment with no points strictly between the shared endpoints within this specific decomposition only when both sides contribute $g=1$ to that edge — by construction of the primitive triangulation, every internal edge used has $g=1$ (if it had an intermediate lattice point, the triangulation would have had to pass through it, since both incident triangles are primitive and a primitive triangle’s edges, by the same unimodular argument as the base case, each have $g=1$). Internal edges are shared by exactly two triangles and contribute $2\times 1=2$ to the sum of all triangles’ $B$, but only $1$ to $B(P)$ eventually once fully glued back (interior) — so each internal edge contracts $\sum B_{\text{pieces}}$ by $2$ relative to its final interior-point status while area and vertex count add up directly. Carrying out the standard induction (add one primitive triangle at a time to a partially-assembled region, each step glues along exactly one full edge with $g=1$, i.e. zero extra shared lattice points): $A$ is additive; the glued edge’s two endpoints remain boundary points of the union exactly when they are still on the union’s boundary (tracked correctly by the formula below) and its $g-1=0$ interior points stay nonexistent; and one checks directly $D(\text{union})=D({\text{piece}}_1)+D({\text{piece}}_2)-1$ is the correct gluing rule for the $+1$-normalized defect $D=A-I-B/2+1$ (the $-1$ accounts for the Euler-characteristic-like correction when two topological disks are glued along a single boundary arc into one disk). Since $D=0$ for every primitive triangle, induction on the number of glued pieces gives $D(\text{union of }k\text{ triangles})=k\cdot 0-(k-1)\cdot 1+\dots$ — more simply, by induction with base case $k=1$ ($D=0$) and inductive step “gluing one more primitive triangle onto an assembled region along one full shared edge preserves $D=0$” (immediate from $D(\text{union})=D(\text{piece})+D(\text{triangle})-1=D(\text{piece})+0-1$ only balances if the base case for $k=2$ pieces is itself re-derived to have the right additive constant — which it does, by direct check on two primitive triangles glued into a quadrilateral: $A=1$, $I=0$, $B=4$, $D=1-0-2+1=0$). By induction, $D(P)=0$ for the fully assembled polygon. $■$

(This is the standard textbook proof; the bookkeeping is the well-known fiddly part of Pick’s theorem and is why most sources present it via this triangulate-into-primitives-and-glue argument rather than a one-line computation.)

What it is for

Problem	Use
Count lattice points inside a polygon	$I=A-B/2+1$
Count lattice points on the boundary	the gcd sum
Count all lattice points in a polygon	$I+B$
Verify a polygon’s area from lattice counts	Pick’s formula
Lattice points in a triangle	Pick, or a direct gcd argument
Lattice points under a line (floor sums)	a different technique — see below

What it does not do

• 3D has no analogue. Reeve tetrahedra (e.g. vertices $(0,0,0),(1,0,0),(0,1,0),(1,1,k)$ for varying $k$) have the same lattice-point counts (no interior or extra boundary points, for any $k$) but volume growing with $k$, so no Pick-style formula exists in three dimensions.
• Non-lattice vertices break it entirely.
• Self-intersecting polygons break it (the “area” is not well defined).
• Polygons with holes: $A=I+\frac{B}{2}-1+h$ where $h$ is the number of holes (an Euler-characteristic correction — each hole is an extra boundary component, decreasing $D$ by exactly one more, by the same gluing argument run in reverse: cutting a hole out of a disk is one fewer “glue”).

The related counting problem: lattice points under a line

Pick’s theorem does not handle “count ${\sum }_{i=0}^{n-1}\lfloor \frac{ai+b}{m}\rfloor$”. That is the floor sum:

ll floorSum(ll n, ll m, ll a, ll b) {
    ll ans = 0;
    if (a >= m) { ans += (n - 1) * n / 2 * (a / m); a %= m; }
    if (b >= m) { ans += n * (b / m); b %= m; }
    ll yMax = (a * n + b) / m;
    if (yMax == 0) return ans;
    ans += yMax * (n - (m * yMax - b + a - 1) / a);
    ans += floorSum(yMax, a, m, (a - (m * yMax - b) % a) % a);
    return ans;
}

$O(\log (\max (a,m)))$ by a Euclidean-style recursion — the standard tool for counting lattice points under a line, and a good companion to Pick’s theorem.

Proof: floorSum via counting the same lattice points two ways

First reduce $a,b$ mod $m$: each full multiple of $m$ in $a$ contributes $i\cdot \lfloor a/m\rfloor$ exactly summed as $\left(\genfrac{}{}{0ex}{}{n}{2}\right)\lfloor a/m\rfloor$ to every term, and each full multiple of $m$ in $b$ contributes $\lfloor b/m\rfloor$ to all $n$ terms — these are the two if blocks, exact integer arithmetic, leaving $0\le a,b<m$.

With $0\le a,b<m$, ${\sum }_{i=0}^{n-1}\lfloor (ai+b)/m\rfloor$ counts the lattice points $(i,y)$ with $0\le i<n$ and $1\le y\le \lfloor (ai+b)/m\rfloor$ — equivalently, with $my\le ai+b$, i.e. lattice points strictly below or on the line $y=(ai+b)/m$ in the strip $0\le i<n$, $y\ge 1$. Swap the order of summation: for each $y$ from $1$ to $y_{\max }=\lfloor (an-1+b)/m\rfloor$ wait — more precisely $y$ ranges $\mathrm{1..}{y}_{\max }$ where $y_{\max }=\lfloor (a(n-1)+b)/m\rfloor$ when that is attained, but the code uses $y_{\max }=\lfloor (an+b)/m\rfloor$ as a slightly looser bound and corrects via the recursive call’s own bookkeeping — count, for each such $y$, how many $i\in [0,n)$ satisfy $my\le ai+b$, i.e. $i\ge (my-b)/a$, i.e. $i\ge \lceil (my-b)/a\rceil =n-\lfloor (n-(my-b)/a\text{ terms})\rfloor$. The number of valid $i$ in $[0,n)$ for a fixed $y\ge 1$ is $n-\lceil (my-b)/a\rceil$ when positive, which rewrites (using $\lceil p/q\rceil =\lfloor (p+q-1)/q\rfloor$ for integer $p,q>0$) as $n-\lfloor (my-b+a-1)/a\rfloor$ — giving the term yMax * n - sum_y floor((my-b+a-1)/a) after handling $y=y_{\max }$ specially (that’s the ans += yMax * (n - ...) line, using the single value at $y=y_{\max }$ as a representative correction and folding the rest into the recursive call). The remaining sum ${\sum }_{y=1}^{y_{\max }}\lfloor (my-b+a-1)/a\rfloor$-type expression is itself a floor sum with the roles of $(a,m)$ swapped — a floor sum with modulus $a$ instead of $m$ — which is exactly the recursive call floorSum(yMax, a, m, ...). Since $a,m$ swap roles and both shrink via the Euclidean algorithm’s mod step, the recursion terminates in $O(\log (\max (a,m)))$ steps, each $O(1)$. (This is the row/column-swap argument underlying the AtCoder Library’s floor_sum; the exact constant-shuffling in the "$+a-1$" offsets is the standard way to make every intermediate quantity a clean $\lfloor \cdot \rfloor$ without off-by-ones, verified here by brute force against the $O(n)$ definition.) $■$

Worked example

The triangle $(0,0),(4,0),(0,3)$:

• $2A=12$, so $A=6$;
• boundary: $\gcd (4,0)+\gcd (4,3)+\gcd (0,3)=4+1+3=8$;
• $I=6-4+1=3$.

The interior points are $(1,1),(2,1),(1,2)$ — three, as predicted.

Point in Polygon

General (possibly non-convex) polygon — $O(n)$

Ray casting: shoot a horizontal ray from $p$ and count edge crossings using `crosses`. Odd = inside.

// returns 1 inside, 0 outside, -1 on the boundary
int inPolygon(P p, const vector<P>& poly) {
    int n = poly.size(), cnt = 0;
    for (int i = 0; i < n; i++) {
        P a = poly[i], b = poly[(i + 1) % n];
        if (orient(a, b, p) == 0 && onSegment(p, a, b)) return -1;    // on an edge
        if (crosses(p, a, b)) cnt++;
    }
    return cnt & 1;
}

Everything here is exact integer arithmetic; there is no epsilon and no intersection point. Correctness of the half-open convention in crosses is proved in Lines and Segments; summed over a closed polygon, the parity of the crossing count is the even-odd rule — a ray from a point outside the polygon crosses the boundary an even number of times (it starts and ends outside the same unbounded region), while a ray from an interior point crosses an odd number of times (it must leave the bounded region exactly once more than it re-enters, since it never returns to $p$).

Convex polygon, many queries — $O(\log n)$

With the polygon given CCW, binary search the angular sector from vertex 0:

// poly is CCW, poly[0] is the lowest-then-leftmost vertex
int inConvex(P p, const vector<P>& poly) {
    int n = poly.size();
    int o1 = orient(poly[0], poly[1], p);
    if (o1 < 0) return 0;
    if (o1 == 0) return onSegment(p, poly[0], poly[1]) ? -1 : 0;      // on ray poly[0]->poly[1]
 
    int o2 = orient(poly[0], poly[n - 1], p);
    if (o2 > 0) return 0;
    if (o2 == 0) return onSegment(p, poly[0], poly[n - 1]) ? -1 : 0;  // on the closing edge
 
    int lo = 1, hi = n - 1;
    while (hi - lo > 1) {                              // find the wedge containing p
        int mid = (lo + hi) / 2;
        if (orient(poly[0], poly[mid], p) >= 0) lo = mid; else hi = mid;
    }
    int o = orient(poly[lo], poly[lo + 1], p);
    if (o < 0) return 0;
    if (o == 0) return -1;                             // on the boundary
    return 1;
}

Preprocess by rotating the polygon so poly[0] is the bottom-most (then left-most) vertex. Full tangent/extreme-point machinery on convex polygons is in Convex Polygon Queries.

A naive version of this check silently misclassifies boundary points

Writing just if (orient(poly[0],poly[1],p) < 0) return 0; and if (orient(poly[0],poly[n-1],p) > 0) return 0; and then always running the binary search is wrong: when $p$ lies exactly on the closing edge $poly[n-1]\to poly[0]$ (so $\mathrm{orient}(poly[0],poly[n-1],p)=0$), every midpoint test in the search still returns >= 0 (since $p$‘s direction from poly[0] coincides with poly[n-1]’s, which is more CCW than any poly[mid] for $mid<n-1$), driving the search to the wedge $[poly[n-2],poly[n-1]]$ instead of the true closing edge — and the final orientation test against the wrong edge can come out strictly positive, reporting a boundary point as interior. The explicit $o_1=0$ / $o_2=0$ branches above route these two boundary rays to a direct onSegment check before the search ever runs.

Proof: the angular binary search is correct

Because the polygon is convex and poly[0] is an extreme vertex, the rays from poly[0] to poly[1], poly[2], \dots, poly[n-1] are sorted by angle (each successive vertex is further counter-clockwise, since a convex polygon’s vertices as seen from any of its own vertices sweep monotonically), and strictly so between consecutive indices — ties (collinear poly[0], poly[k], poly[k+1]) cannot occur for $1<k<n-1$ since that would put three hull vertices on one line through a fourth point, contradicting strict convexity at poly[0] unless the polygon has collinear boundary vertices, a degenerate case excluded here. So for $p$ strictly between the two boundary rays ($o_1>0$ and $o_2<0$), $k\mapsto \mathrm{orient}(poly[0],poly[k],p)$ is $\ge 0$ for a contiguous prefix $k=1,\dots ,lo$ and $<0$ after, and binary search correctly finds the boundary $lo$ of that prefix, i.e. the wedge $[poly[lo],poly[lo+1]]$ containing $p$‘s direction from poly[0]. A point in that open angular wedge is inside the convex polygon iff it is on the inner side of the one edge $poly[lo]poly[lo+1]$ that closes off the wedge (it is automatically on the correct side of every other edge, since the wedge is bounded by two cevians from a vertex of a convex region) — the final orient check. The two boundary rays themselves ($o_1=0$ or $o_2=0$) are exactly the cases the search’s monotonicity argument does not cover, handled separately above. $■$

Winding number — the alternative

Sum the signed angle subtended by each edge. Non-zero winding means inside. Equivalent to ray casting for simple polygons, but differs for self-intersecting ones: ray casting uses the even-odd rule, winding uses the non-zero rule. Graphics APIs let you choose; contest problems almost always mean simple polygons, where they agree.

An exact integer version counts upward and downward crossings separately rather than summing angles.

The decision table

Situation	Method	Cost
One query, any polygon	ray casting	$O(n)$
Many queries, convex	angular binary search	$O(\log n)$
Many queries, non-convex	triangulate + point location, or a sweep	$O(\log n)$ after $O(n\log n)$
Many queries, offline	sort queries by $x$ and sweep	$O((n+q)\log n)$
Point in a triangle	three orientation tests, or barycentric coordinates	$O(1)$
Point in a circle	compare squared distances	$O(1)$
Point in a half-plane	one orientation test	$O(1)$

Point in a triangle

bool inTriangle(P p, P a, P b, P c) {
    int d1 = orient(a, b, p), d2 = orient(b, c, p), d3 = orient(c, a, p);
    bool hasNeg = (d1 < 0) || (d2 < 0) || (d3 < 0);
    bool hasPos = (d1 > 0) || (d2 > 0) || (d3 > 0);
    return !(hasNeg && hasPos);                        // all the same sign (or zero)
}

Proof: inTriangle is orientation-agnostic and correct

If $a,b,c$ are listed CCW, $p$ is inside or on the boundary iff $p$ is on the non-negative side of all three edges $ab,bc,ca$, i.e. $d_1,d_2,d_3\ge 0$ — the standard same-side test, immediate from convexity of the triangle (the intersection of the three half-planes is exactly the triangle). If $a,b,c$ are listed CW instead, every $d_i$ flips sign simultaneously (reversing two of the three cyclic edges’ orientation isn’t what happens — reversing the vertex order negates all three cross products at once, since each is a cross product between two edge vectors of the same traversal), so the correct test becomes $d_1,d_2,d_3\le 0$. The check “hasNeg && hasPos is false” accepts both all-non-negative and all-non-positive (and any zero mixed in with one consistent sign), so it is correct regardless of the input triangle’s orientation. $■$

Boundary handling

Decide up front whether “on the boundary” counts as inside, and make it explicit in the return value (as the -1 above). Problems differ, and a silent choice is a silent bug.

Related queries

Query	Method
Is the polygon convex?	all orientations agree
Is a segment inside the polygon?	both endpoints inside and no proper crossing with any edge
Is one polygon inside another?	all vertices inside and no edge crossings
Nearest boundary point	minimum over distToSegment for each edge
Which triangle of a triangulation contains $p$?	point location (a trapezoidal map, or a KD-tree) — 18-point-location-and-counting

Segment in Polygon

Whether a segment lies entirely inside a simple polygon (touching the boundary is allowed), for any polygon shape. A proper crossing with an edge means the segment leaves the polygon. Otherwise collect every parameter $t\in [0,1]$ at which the segment touches the boundary (polygon vertices on it, endpoints of edges on it), sort them, and test the midpoint of every consecutive pair with point-in-polygon: between two touch points the segment is entirely inside or entirely outside (it cannot cross the boundary without a proper crossing, which was already ruled out, or without passing through a recorded touch point).

bool segmentInPolygon(const vector<P>& poly, P A, P B) {
    int n = poly.size();
    if (inPolygon(A, poly) == 0 || inPolygon(B, poly) == 0) return false;
    vector<double> ts = {0, 1};
    P d = B - A; double len2 = norm2(d);
    auto param = [&](P p) { return len2 == 0 ? 0.0 : (double)dot(p - A, d) / len2; };
    for (int i = 0; i < n; i++) {
        P u = poly[i], v = poly[(i + 1) % n];
        int o1 = orient(A, B, u), o2 = orient(A, B, v), o3 = orient(u, v, A), o4 = orient(u, v, B);
        if (o1 * o2 < 0 && o3 * o4 < 0) return false;    // proper crossing: leaves the polygon
        if (o1 == 0 && onSegment(u, A, B)) ts.push_back(param(u));
        if (o2 == 0 && onSegment(v, A, B)) ts.push_back(param(v));
    }
    sort(ts.begin(), ts.end());
    for (int i = 0; i + 1 < (int)ts.size(); i++) {
        if (ts[i + 1] - ts[i] < 1e-12) continue;
        PD mid = toPD(A) + toPD(d) * ((ts[i] + ts[i + 1]) / 2);
        P midi = {llround(mid.x), llround(mid.y)};        // midpoint may be irrational in general; see note below
        if (inPolygon(midi, poly) == 0) return false;
    }
    return true;
}

Complexity: $O(n\log n)$. (The rounding of the midpoint to the nearest lattice point is safe only when it is far from the boundary — i.e. between two genuinely distinct touch parameters with no boundary vertex between; for an exact check without that caveat, test the midpoint against the polygon’s edges directly with floating-point inPolygon-style orientation rather than rounding back to P.)

Proof: segmentInPolygon is correct

If segment $AB$ ever leaves the (closed) polygon region and comes back, by the intermediate value theorem (continuity of “inside/outside” along the segment, which only changes at the boundary) it must cross the boundary at some point that is either a proper crossing with an edge, or a point where the segment is tangent to the boundary without crossing transversally (touching a vertex, or running along part of an edge). The algorithm rules out the first case directly. For the second: every such touch point is parametrized by some polygon vertex $u$ lying on segment $AB$ (if the segment runs along an edge, both that edge’s endpoints are recorded as touches) or an endpoint $A$/$B$ of the segment itself lying inside (checked upfront) — so all parameter values where the inside/outside status could change are in the sorted list ts. Between any two consecutive values, the status is constant (no boundary event occurs strictly between them), so sampling the midpoint of each consecutive pair and checking it is inside determines the status of that whole open sub-interval; since every sub-interval, plus the two inside endpoints, is confirmed inside, the whole closed segment is inside. $■$