7.1.2 Line

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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179

/*

A straight line in two dimensions, templated on the coordinate type `T` (default `double`). The line
$ax + by + c = 0$ is stored unnormalized, so a line through integer points $(p_x, p_y)$ and
$(q_x, q_y)$ keeps the integer coefficients: $a = q_y - p_y$, $b = p_x - q_x$,
$c = -(a \cdot p_x + b \cdot p_y)$. Storage and the exact predicates `contains()`, `is_parallel()`,
and `is_perpendicular()` therefore stay exact for integral type `T`. The inherently-fractional
operations `slope()`, `x()`, and `y()` return `fp_t`, which is `double` when `T` is integral and `T`
itself otherwise.

Because coefficients are not normalized to a canonical form, `operator==` compares lines by
proportionality (whether they are the same geometric line) rather than by coefficient equality.

Operations include `horizontal()`, `vertical()`, `slope()`, evaluating $y$ at some $x$ with `y()`
(and vice versa with `x()`), checking if a point falls on the line with `contains()`, checking if
another line is parallel or perpendicular with `is_parallel()` or `is_perpendicular()`, and finding
the `parallel()` or `perpendicular()` line through a point.

Overflow warning: the exact predicates form products of coefficients (e.g. cross terms in
`is_parallel()`/`is_perpendicular()` and `a*p.x + b*p.y` in `contains()`), which grow like the
squared coordinate magnitude. For integral type `T`, use `LineL` (`line<int64_t>`) once
coordinates exceed ~46000, or the 32-bit products overflow.

Type aliases:
- `LineI = line<int>`: exact integer-coefficient lines (small values only; see overflow warning)
- `LineL = line<int64_t>`: exact integer-coefficient lines for large coordinates
- `LineD = line<double>`: standard floating-point
- `LineLD = line<long double>`: extra precision
- `Line = LineD`: default line type is double

Time Complexity:
- O(1) per call to the constructor and all other operations.

Space Complexity:
- O(1) for storage of the line.
- O(1) auxiliary for all operations.

*/

#include <cmath>
#include <cstdint>
#include <limits>
#include <ostream>
#include <type_traits>
#include <utility>

const double EPS = 1e-9;
const double M_NAN = std::numeric_limits<double>::quiet_NaN();

// SFINAE guard: valid only when Pt exposes numeric .x/.y members, so templated point constructors
// don't hijack scalar-argument calls.
template<class Pt>
using if_point = decltype(std::declval<const Pt &>().x, std::declval<const Pt &>().y, void());

template<class T>
struct line {
  // fp_t is T for floating-point types, double for integral types.
  using fp_t = std::conditional_t<std::is_integral<T>::value, double, T>;

  T a, b, c;

  line() : a(0), b(0), c(0) {}  // Invalid or uninitialized line.
  line(T a, T b, T c) : a(a), b(b), c(c) {}

  // Line through two points. Coefficients are exact (type T) for integer points.
  template<class Pt, class = if_point<Pt>>
  line(const Pt &p, const Pt &q) : a(q.y - p.y), b(p.x - q.x), c(-(a * p.x + b * p.y)) {}

  // Line of a given slope through a point. Intended for floating-point T (a generic slope
  // produces non-integer coefficients), so this constructor is disabled for integer T.
  template<
      class Pt, class = if_point<Pt>, class U = T,
      class = std::enable_if_t<std::is_floating_point<U>::value>>
  line(fp_t slope, const Pt &p) : a(-slope), b(1), c(slope * p.x - p.y) {}

  // Epsilon-aware for floating-point T; exact for int and for coordinate types like Modular or
  // Rational, which therefore compose for all of the predicates below.
  template<class A, class B>
  static bool eq(const A &p, const B &q) {
    if constexpr (std::is_floating_point<T>::value) {
      return fabs(p - q) <= EPS;
    }
    return p == q;
  }

  // Two lines are equal iff their coefficient vectors are proportional (the same line).
  bool operator==(const line &l) const {
    return eq(a * l.b, l.a * b) && eq(b * l.c, l.b * c) && eq(a * l.c, l.a * c);
  }

  bool operator!=(const line &l) const { return !(*this == l); }

  // Returns whether the line is a valid line (not both a and b zero).
  bool valid() const { return !(eq(a, 0) && eq(b, 0)); }
  bool horizontal() const { return valid() && eq(a, 0); }
  bool vertical() const { return valid() && eq(b, 0); }

  // Slope -a/b, or NaN if the line is vertical or invalid.
  fp_t slope() const { return (!valid() || eq(b, 0)) ? M_NAN : -(fp_t)a / b; }

  // Solve for x at a given y (NaN if horizontal or invalid).
  fp_t x(fp_t y) const { return (!valid() || eq(a, 0)) ? M_NAN : -((fp_t)b * y + c) / a; }

  // Solve for y at a given x (NaN if vertical or invalid).
  fp_t y(fp_t x) const { return (!valid() || eq(b, 0)) ? M_NAN : -((fp_t)a * x + c) / b; }

  // Whether point p lies on the line. Exact for integer T and integer p.
  template<class Pt>
  bool contains(const Pt &p) const {
    return eq(a * p.x + b * p.y + c, 0);
  }

  // Parallel iff the normals are parallel (cross == 0); perpendicular iff normals are
  // perpendicular (dot == 0). Both exact for integer T.
  bool is_parallel(const line &l) const { return eq(a * l.b, l.a * b); }
  bool is_perpendicular(const line &l) const { return eq(a * l.a, -(b * l.b)); }

  // Parallel/perpendicular line through point p. Exact for integer T and integer p.
  template<class Pt, class = if_point<Pt>>
  line parallel(const Pt &p) const {
    return line(a, b, -(a * p.x + b * p.y));
  }

  template<class Pt, class = if_point<Pt>>
  line perpendicular(const Pt &p) const {
    return line(-b, a, b * p.x - a * p.y);
  }

  friend std::ostream &operator<<(std::ostream &out, const line &l) {
    auto clean = [](T v) -> double {
      return fabs(static_cast<double>(v)) < EPS ? 0 : static_cast<double>(v);
    };
    return out << clean(l.a) << "x" << std::showpos << clean(l.b) << "y" << clean(l.c) << "=0"
               << std::noshowpos;
  }
};

using LineI = line<int>;
using LineL = line<int64_t>;
using LineD = line<double>;
using LineLD = line<long double>;
using Line = LineD;  // Default line type is double.

/*** Example Usage ***/

#include <cassert>

struct Point {
  double x, y;
  Point(double x = 0, double y = 0) : x(x), y(y) {}
};

struct PointI {
  int x, y;
  PointI(int x = 0, int y = 0) : x(x), y(y) {}
};

int main() {
  // Floating-point line.
  Line l(2, -5, -8);
  Line para = Line(2, -5, -8).parallel(Point(-6, -2));
  Line perp = Line(2, -5, -8).perpendicular(Point(-6, -2));
  assert(l.is_parallel(para) && l.is_perpendicular(perp));
  assert(l.slope() == 0.4);
  // operator== is proportionality-based: para is the same line as -0.4x + y - 0.4 = 0.
  assert(para == Line(-0.4, 1, -0.4));
  assert(perp == Line(2.5, 1, 17));
  assert(l.contains(Point(1.5, -1)));  // 2(1.5) - 5(-1) - 8 = 0.

  // Integer line through integer points - exact coefficients and predicates.
  LineI il(PointI(0, 0), PointI(2, 1));  // a=1, b=-2, c=0  ->  x - 2y = 0
  assert(il.a == 1 && il.b == -2 && il.c == 0);
  assert(il.contains(PointI(4, 2)));  // exact integer check
  assert(!il.contains(PointI(4, 3)));
  assert(il.is_parallel(LineI(PointI(1, 1), PointI(3, 2))));        // both slope 1/2
  assert(il.is_perpendicular(LineI(PointI(0, 0), PointI(1, -2))));  // dot of normals = 0
  assert(il.slope() == 0.5);                                        // fp_t result
  return 0;
}