Alex's Anthology of Algorithms Common Code for Contests in Concise C++
Geometry / Geometry Primitives

A 2D point class templated on the coordinate type T. Algebraic operations preserve T, while metric operations use fp_t, which is T for floating-point coordinates and double otherwise.

Non-floating-point coordinate types such as Modular or Rational compose for exact operations (arithmetic, dot, cross, sqnorm, comparisons, and EQ) using the coordinate type's own operators. Their metric operations convert coordinates explicitly to double; for Modular, this uses the stored representative and is not finite-field geometry.

  • TPoint<T>() constructs the origin, while TPoint<T>(x, y) constructs point $(x, y)$. A pair of values can also be converted explicitly to a point.
  • Operators +, -, *, /, and their compound forms act element-wise on points or apply a scalar. Comparisons are exact and lexicographic, as required by standard algorithms and containers. EQ(p, q) instead uses EPS for floating-point coordinates and remains exact for other coordinate types.
  • sqnorm(), dot(p), and cross(p) return the exact algebraic results in coordinate type T.
  • norm(), arg(), proj(p), and normalize() return metric results in fp_t; division also promotes to fp_t to avoid integer truncation.
  • rotate90(), rotate180(), and rotate270() rotate by the corresponding counter-clockwise cardinal angle; overloads taking p rotate about point p.
  • rotate_cw(t) and rotate_ccw(t) rotate by an arbitrary angle t in radians; overloads taking (p, t) rotate about point p.
  • reflect(p) reflects across point p, while reflect(p, q) reflects across the line through points p and q.
  • to_double() and to_ldouble() explicitly convert the coordinate type. Integral points also convert implicitly to floating-point point types, so PointI can be passed where PointD is expected.
  • PointI, PointL, PointD, and PointLD use int, long long, double, and long double coordinates, respectively. Point aliases PointD.

Overflow warning: the exact products dot(), cross(), and sqnorm() grow like the squared coordinate magnitude. With PointI these overflow a 32-bit int once coordinates exceed roughly a few tens of thousands, so use PointL for larger integer coordinates.

Implementation

#include <cmath>
#include <ostream>
#include <tuple>
#include <type_traits>
#include <utility>

const double EPS = 1e-9;

// Epsilon-aware for floating-point coordinates; exact for int and for coordinate types like Modular
// or Rational, which therefore compose for all of the predicates below.
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);
}

template<typename T>
struct TPoint {
  // Metric operations preserve floating-point coordinates and otherwise promote to double.
  using fp_t = std::conditional_t<std::is_floating_point<T>::value, T, double>;

  T x, y;

  TPoint() : x(0), y(0) {}
  TPoint(T x, T y) : x(x), y(y) {}
  explicit TPoint(const std::pair<T, T> &p) : x(p.first), y(p.second) {}

  // Implicit conversion from integral to floating-point point (optional, can just use to_double()).
  template<
      typename U,
      typename = std::enable_if_t<std::is_integral<U>::value && std::is_floating_point<T>::value>>
  TPoint(const TPoint<U> &p) : x(static_cast<T>(p.x)), y(static_cast<T>(p.y)) {}

  friend bool EQ(const TPoint &a, const TPoint &b) { return EQ(a.x, b.x) && EQ(a.y, b.y); }

  bool operator==(const TPoint &p) const { return x == p.x && y == p.y; }
  bool operator!=(const TPoint &p) const { return !(*this == p); }
  bool operator<(const TPoint &p) const { return std::tie(x, y) < std::tie(p.x, p.y); }
  bool operator>(const TPoint &p) const { return p < *this; }
  bool operator<=(const TPoint &p) const { return !(*this > p); }
  bool operator>=(const TPoint &p) const { return !(*this < p); }
  TPoint operator+(const TPoint &p) const { return {x + p.x, y + p.y}; }
  TPoint operator-(const TPoint &p) const { return {x - p.x, y - p.y}; }
  TPoint operator+(T v) const { return {x + v, y + v}; }
  TPoint operator-(T v) const { return {x - v, y - v}; }
  TPoint operator*(T v) const { return {x * v, y * v}; }

  // Division always promotes to fp_t to avoid integer truncation.
  TPoint<fp_t> operator/(fp_t v) const { return {fp_t(x) / v, fp_t(y) / v}; }

  TPoint &operator+=(const TPoint &p) { x += p.x; y += p.y; return *this; }
  TPoint &operator-=(const TPoint &p) { x -= p.x; y -= p.y; return *this; }
  TPoint &operator+=(T v) { x += v; y += v; return *this; }
  TPoint &operator-=(T v) { x -= v; y -= v; return *this; }
  TPoint &operator*=(T v) { x *= v; y *= v; return *this; }
  friend TPoint operator+(T v, const TPoint &p) { return p + v; }
  friend TPoint operator*(T v, const TPoint &p) { return p * v; }

  // --- Exact operations: return T or TPoint<T>, work for any coordinate type ---

  T sqnorm() const { return x * x + y * y; }                    // Overflow warning.
  T dot(const TPoint &p) const { return x * p.x + y * p.y; }    // Overflow warning.
  T cross(const TPoint &p) const { return x * p.y - y * p.x; }  // Overflow warning.

  // --- Floating-point operations: return fp_t or TPoint<fp_t> ---

  fp_t norm() const { return std::hypot(fp_t(x), fp_t(y)); }
  fp_t arg() const { return std::atan2(fp_t(y), fp_t(x)); }
  fp_t proj(const TPoint &p) const { return (fp_t(x) * p.x + fp_t(y) * p.y) / p.norm(); }

  TPoint<fp_t> normalize() const {
    fp_t n = norm();
    if (n < 1e-30) {              // guard against dividing by a near-zero norm, not a
      return {fp_t{0}, fp_t{0}};  // geometric tolerance (EPS is too coarse here)
    }
    return {fp_t(x) / n, fp_t(y) / n};
  }

  // --- Cardinal rotations: exact for all coordinate types including int ---

  // Returns (x, y) rotated 90/180/270 degrees counter-clockwise about the origin.
  TPoint rotate90() const { return {-y, x}; }
  TPoint rotate180() const { return {-x, -y}; }
  TPoint rotate270() const { return {y, -x}; }

  // Returns (x, y) rotated 90/180/270 degrees counter-clockwise about point p.
  TPoint rotate90(const TPoint &p) const { return (*this - p).rotate90() + p; }
  TPoint rotate180(const TPoint &p) const { return (*this - p).rotate180() + p; }
  TPoint rotate270(const TPoint &p) const { return (*this - p).rotate270() + p; }

  // --- Arbitrary-angle rotations: always return floating-point ---

  // Returns (x, y) rotated t radians clockwise about the origin.
  TPoint<fp_t> rotate_cw(fp_t t) const {
    fp_t fx(x), fy(y);
    return {fx * std::cos(t) + fy * std::sin(t), fy * std::cos(t) - fx * std::sin(t)};
  }

  // Returns (x, y) rotated t radians counter-clockwise about the origin.
  TPoint<fp_t> rotate_ccw(fp_t t) const {
    fp_t fx(x), fy(y);
    return {fx * std::cos(t) - fy * std::sin(t), fx * std::sin(t) + fy * std::cos(t)};
  }

  // Returns (x, y) rotated t radians clockwise about point p.
  TPoint<fp_t> rotate_cw(const TPoint &p, fp_t t) const {
    return TPoint<fp_t>{fp_t(x) - p.x, fp_t(y) - p.y}.rotate_cw(t) +
           TPoint<fp_t>{fp_t(p.x), fp_t(p.y)};
  }

  // Returns (x, y) rotated t radians counter-clockwise about point p.
  TPoint<fp_t> rotate_ccw(const TPoint &p, fp_t t) const {
    return TPoint<fp_t>{fp_t(x) - p.x, fp_t(y) - p.y}.rotate_ccw(t) +
           TPoint<fp_t>{fp_t(p.x), fp_t(p.y)};
  }

  // --- Reflections ---

  // Returns (x, y) reflected across point p. Exact for any coordinate type.
  TPoint reflect(const TPoint &p) const { return {2 * p.x - x, 2 * p.y - y}; }

  // Returns (x, y) reflected across the line containing points p and q.
  // Always returns floating-point coordinates.
  TPoint<fp_t> reflect(const TPoint &p, const TPoint &q) const {
    TPoint<fp_t> fp{fp_t(p.x), fp_t(p.y)};
    if (EQ(p, q)) {
      return TPoint<fp_t>{fp_t(x), fp_t(y)}.reflect(fp);
    }
    TPoint<fp_t> r{fp_t(x) - p.x, fp_t(y) - p.y};
    TPoint<fp_t> s{fp_t(q.x) - p.x, fp_t(q.y) - p.y};
    fp_t ssq = s.sqnorm();
    r = TPoint<fp_t>{(r.x * s.x + r.y * s.y) / ssq, (r.x * s.y - r.y * s.x) / ssq};
    return TPoint<fp_t>{r.x * s.x - r.y * s.y + fp.x, r.x * s.y + r.y * s.x + fp.y};
  }

  // --- Explicit type conversions ---

  TPoint<double> to_double() const { return {static_cast<double>(x), static_cast<double>(y)}; }

  TPoint<long double> to_ldouble() const {
    return {static_cast<long double>(x), static_cast<long double>(y)};
  }

  // --- Friend free-function versions ---

  friend T sqnorm(const TPoint &p) { return p.sqnorm(); }
  friend fp_t norm(const TPoint &p) { return p.norm(); }
  friend fp_t arg(const TPoint &p) { return p.arg(); }
  friend T dot(const TPoint &p, const TPoint &q) { return p.dot(q); }
  friend T cross(const TPoint &p, const TPoint &q) { return p.cross(q); }
  friend fp_t proj(const TPoint &p, const TPoint &q) { return p.proj(q); }
  friend TPoint<fp_t> normalize(const TPoint &p) { return p.normalize(); }
  friend TPoint rotate90(const TPoint &p) { return p.rotate90(); }
  friend TPoint rotate180(const TPoint &p) { return p.rotate180(); }
  friend TPoint rotate270(const TPoint &p) { return p.rotate270(); }
  friend TPoint rotate90(const TPoint &p, const TPoint &q) { return p.rotate90(q); }
  friend TPoint rotate180(const TPoint &p, const TPoint &q) { return p.rotate180(q); }
  friend TPoint rotate270(const TPoint &p, const TPoint &q) { return p.rotate270(q); }
  friend TPoint<fp_t> rotate_cw(const TPoint &p, fp_t t) { return p.rotate_cw(t); }
  friend TPoint<fp_t> rotate_ccw(const TPoint &p, fp_t t) { return p.rotate_ccw(t); }
  friend TPoint<fp_t> rotate_cw(const TPoint &p, const TPoint &q, fp_t t) { return p.rotate_cw(q, t); }
  friend TPoint<fp_t> rotate_ccw(const TPoint &p, const TPoint &q, fp_t t) { return p.rotate_ccw(q, t); }
  friend TPoint reflect(const TPoint &p, const TPoint &q) { return p.reflect(q); }
  friend TPoint<fp_t> reflect(const TPoint &p, const TPoint &a, const TPoint &b) { return p.reflect(a, b); }

  friend std::ostream &operator<<(std::ostream &out, const TPoint &p) {
    if constexpr (std::is_floating_point<T>::value) {
      return out << "(" << (std::fabs(p.x) < EPS ? 0 : p.x) << ","
                 << (std::fabs(p.y) < EPS ? 0 : p.y) << ")";
    }
    return out << "(" << p.x << "," << p.y << ")";
  }
};

using PointI = TPoint<int>;
using PointL = TPoint<long long>;
using PointD = TPoint<double>;
using PointLD = TPoint<long double>;
using Point = PointD;  // Default point type is double.

// Can compose with numerical types from chapter 6:
// using PointB = TPoint<BigInt>;
// using PointR = TPoint<Rational<int64_t>>;
// using PointM = TPoint<Modular<1000000007>>;

Example Usage

#include <cassert>
using namespace std;

const double PI = acos(-1.0);

int main() {
  Point p(-10, 3);
  assert(EQ(Point(-18, 29), p + Point(-3, 9) * 6.0 / 2.0 - Point(-1, 1)));
  assert(EQ(109, p.sqnorm()));
  assert(EQ(10.44030650891, p.norm()));
  assert(EQ(2.850135859112, p.arg()));
  assert(EQ(0, p.dot(Point(3, 10))));
  assert(EQ(0, p.cross(Point(10, -3))));
  assert(EQ(10, p.proj(Point(-10, 0))));
  assert(EQ(1, p.normalize().norm()));
  assert(EQ(Point(-3, -10), p.rotate90()));
  assert(EQ(Point(10, -3), p.rotate180()));
  assert(EQ(Point(3, 10), p.rotate270()));
  assert(EQ(Point(3, 12), p.rotate_cw(Point(1, 1), PI / 2)));
  assert(EQ(Point(1, -10), p.rotate_ccw(Point(2, 2), PI / 2)));
  assert(EQ(Point(10, -3), p.reflect(Point(0, 0))));
  assert(EQ(Point(-10, -3), p.reflect(Point(-2, 0), Point(5, 0))));

  // Integer point - exact arithmetic, float-only ops return PointD.
  PointI a(3, 4), b(1, 0);
  assert(a + b == PointI(4, 4));
  assert(a - b == PointI(2, 4));
  assert(a * 2 == PointI(6, 8));
  assert(a.sqnorm() == 25);
  assert(a.dot(b) == 3);
  assert(a.cross(b) == -4);
  assert(EQ(a.norm(), 5.0));
  assert(a.rotate90() == PointI(-4, 3));
  assert(a.rotate180() == PointI(-3, -4));
  assert(a.rotate270() == PointI(4, -3));
  PointD anorm = a.normalize();  // returns PointD
  assert(EQ(anorm.norm(), 1.0));
  PointD adiv = a / 2.0;  // returns PointD (division always promotes)
  assert(EQ(adiv.x, 1.5) && EQ(adiv.y, 2.0));
  auto rotated = PointI(1, 0).rotate_ccw(PI / 4);  // returns PointD (arbitrary rotation promotes)
  double root2over2 = sqrt(2.0) / 2.0;
  assert(EQ(rotated.x, root2over2) && EQ(rotated.y, root2over2));

  // Implicit conversion PointI -> PointD.
  PointD promoted = PointI(1, 2);
  assert(promoted == PointD(1, 2));
  return 0;
}