7.2.1 Angles

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153

/*

Angle calculations in two dimensions. All operations are inherently floating-point (`atan2`, `acos`,
trigonometry), so these functions use the local double-coordinate `PointD` type. Convert integer
points to `PointD` when asking for angles. The constants `DEG` and `RAD` may be used as multipliers
to convert between degrees and radians. For example, if `t` is a value in degrees, then `t * DEG` is
the equivalent angle in radians; if `t` is in radians, then `t * RAD` is the equivalent angle in
degrees.

- `reduce_deg(t)` takes an angle `t` degrees and returns an equivalent angle in the range $[0, 360)$
  degrees (e.g. $-630$ becomes $90$).
- `reduce_rad(t)` takes an angle `t` radians and returns an equivalent angle in the range
  $[0, 2\pi)$ radians (e.g. $720.5$ becomes $0.5$).
- `polar_point(r, t)` returns a two-dimensional Cartesian point given radius `r` and angle `t`
  radians in polar coordinates (see `std::polar()`).
- `polar_angle(p)` returns the angle in radians of the line segment from $(0, 0)$ to point `p`,
  relative counterclockwise to the positive $x$-axis.
- `angle(a, o, b)` returns the smallest angle in radians formed by the points `a`, `o`, `b` with
  vertex at point `o`.
- `angle_between(a, b)` returns the directed counter-clockwise angle in $[0, 2\pi)$ from vector `a`
  to vector `b`, treating both points as vectors from the origin.
- `angle_between(a1, b1, a2, b2)` returns the smaller angle in radians between two lines
  $a_1 x + b_1 y + c_1 = 0$ and $a_2 x + b_2 y + c_2 = 0$, limited to $[0, \pi / 2]$.
- `cross(a, b, o = Pt(0, 0))` returns the signed $z$-component of the three-dimensional cross
  product between points `a` and `b`, where the $z$-coordinates are implicitly zero and the origin
  is shifted to point `o`. This is also double the signed area of the triangle from these three
  points.
- `turn(a, o, b)` returns $1$ if the path `a` $\to$ `o` $\to$ `b` forms a left turn on the plane,
  $0$ if the path forms a straight line segment, or $-1$ if it forms a right turn.

The points passed to `angle()` must differ from its vertex `o`, vectors passed to `angle_between()`
must be nonzero, and coefficient pairs passed to its line overload must represent valid lines.

Overflow warning: `cross()` and `turn()` preserve the point's coordinate type, so their products may
overflow for integral points. Use a 64-bit coordinate type for large integer coordinates.

Time Complexity:
- O(1) per operation.

Space Complexity:
- O(1) auxiliary for all operations.

*/

#include <algorithm>
#include <cmath>
#include <type_traits>

const double EPS = 1e-9;
const double PI = acos(-1.0);
const double DEG = PI / 180;
const double RAD = 180 / PI;

template<typename T, typename U, typename C = std::common_type_t<T, U>>
bool EQ(T a, U b) {
  if constexpr (std::is_floating_point_v<C>) return C(a) == C(b) || std::fabs(C(a) - C(b)) <= EPS;
  return C(a) == C(b);
}

template<typename T, typename U, typename C = std::common_type_t<T, U>>
bool LT(T a, U b) {
  if constexpr (std::is_floating_point_v<C>) return C(a) < C(b) - EPS;
  return C(a) < C(b);
}

double reduce_deg(double t) {
  if (t < -360) {
    return reduce_deg(fmod(t, 360));
  }
  if (t < 0) {
    return t + 360;
  }
  return (t >= 360) ? fmod(t, 360) : t;
}

double reduce_rad(double t) {
  if (t < -2 * PI) {
    return reduce_rad(fmod(t, 2 * PI));
  }
  if (t < 0) {
    return t + 2 * PI;
  }
  return (t >= 2 * PI) ? fmod(t, 2 * PI) : t;
}

template<typename OutPt>
OutPt polar_point(double r, double t) {
  return OutPt(r * cos(t), r * sin(t));
}

template<typename Pt>
double polar_angle(const Pt &p) {
  double t = atan2(static_cast<double>(p.y), static_cast<double>(p.x));
  return (t < 0) ? (t + 2 * PI) : t;
}

template<typename Pt>
double angle(const Pt &a, const Pt &o, const Pt &b) {
  double ux = o.x - a.x, uy = o.y - a.y;
  double vx = o.x - b.x, vy = o.y - b.y;
  double cosine = (ux * vx + uy * vy) / (std::hypot(ux, uy) * std::hypot(vx, vy));
  return acos(std::max(-1.0, std::min(1.0, cosine)));
}

template<typename Pt>
double angle_between(const Pt &a, const Pt &b) {
  double t = atan2(a.x * b.y - a.y * b.x, a.x * b.x + a.y * b.y);
  return (t < 0) ? (t + 2 * PI) : t;
}

double angle_between(double a1, double b1, double a2, double b2) {
  double t = atan2(a1 * b2 - a2 * b1, a1 * a2 + b1 * b2);
  if (t < 0) {
    t += PI;
  }
  return LT(PI / 2, t) ? (PI - t) : t;
}

template<typename Pt>
auto cross(const Pt &a, const Pt &b, const Pt &o = Pt(0, 0)) {
  return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

template<typename Pt>
int turn(const Pt &a, const Pt &o, const Pt &b) {
  auto c = cross(a, b, o);
  return LT(c, 0) ? 1 : (LT(0, c) ? -1 : 0);
}

/*** Example Usage ***/

struct Point {
  double x, y;
  Point(double x = 0, double y = 0) : x(x), y(y) {}
  bool operator==(const Point &p) const { return x == p.x && y == p.y; }
  friend bool EQ(const Point &a, const Point &b) { return EQ(a.x, b.x) && EQ(a.y, b.y); }
};

#include <cassert>

int main() {
  assert(EQ(123, reduce_deg(-8 * 360 + 123)));
  assert(EQ(1.2345, reduce_rad(2 * PI * 8 + 1.2345)));
  assert(EQ(polar_point<Point>(4, PI), Point(-4, 0)));
  assert(EQ(polar_point<Point>(4, -PI / 2), Point(0, -4)));
  assert(EQ(45, polar_angle(Point(5, 5)) * RAD));
  assert(EQ(135 * DEG, polar_angle(Point(-4, 4))));
  assert(EQ(90 * DEG, angle(Point(5, 0), Point(0, 5), Point(-5, 0))));
  assert(EQ(225 * DEG, angle_between(Point(0, 5), Point(5, -5))));
  assert(-1 == cross(Point(0, 1), Point(1, 0), Point(0, 0)));
  assert(-1 == turn(Point(0, 1), Point(0, 0), Point(-5, -5)));
  return 0;
}