lib_defs.hpp

#ifndef E_SMART_LIB_DEFS
#define E_SMART_LIB_DEFS

#include <algorithm>
#include <complex>
#include <cstdio>
#include <initializer_list>
#include <iostream>
#include <iterator>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <stack>
#include <vector>

#define REPI(i, n) for (size_t i = 0; i < (n); ++i)
#define REPD(i, n) for (size_t i = (n); i >= 0; --i)
#define FOR(i, c) for (auto i : (c))
#define ALL(c) (c).begin(), (c).end()
#define LEN(c) ((c).size())

#define IO_SETUP                                                               \
  std::cin.tie(nullptr);                                                       \
  std::ios::sync_with_stdio(false);

#define SEE(i)                                                                 \
  std::cerr << #i << ": " << i << " at: (" << __FILE__ << ":" << __LINE__      \
            << ")\n";

template <class T, class Comp>
inline void more_impl(T const &src, T &target, Comp f) {
  if (f(target, src)) {
    target = src;
  }
}

template <class T> inline void more_greater(T const &src, T &target) {
  more_impl<T>(src, target, std::greater<T>());
}

template <class T> inline void more_less(T const &src, T &target) {
  more_impl<T>(src, target, std::less<T>());
}

template <class T> inline T diff_abs(T const &l, T const &r) {
  return std::max(l, r) - std::min(l, r);
}

#endif // E_SMART_LIB_DEFS

lib_geometry.hpp

#ifndef E_SMART_LIB_GEOMETRY
#define E_SMART_LIB_GEOMETRY

#include "lib_defs.hpp"

struct point;
struct line;
struct segment;
struct circle;

struct point : public std::complex<fuzzy> {
  point() : std::complex<fuzzy>(0, 0) {}
  explicit point(fuzzy r, fuzzy i) : std::complex<fuzzy>(r, i) {}
  point(std::complex<fuzzy> const &a) : std::complex<fuzzy>(a) {}

  point(point const &old) : std::complex<fuzzy>(old.real(), old.imag()) {}

  fuzzy dot(point const &b) const { return (b * point(conj(*this))).real(); }
  fuzzy cross(point const &b) const { return (b * point(conj(*this))).imag(); }

  int ccw(point const &b, point const &c) const {
    point x = b - *this;
    point y = c - *this;
    if (x.cross(y) > 0)
      return +1;
    if (x.cross(y) < -0)
      return -1;
    if (x.dot(y) < -0)
      return +2;
    if (norm(x) < norm(y))
      return -2;
    return 0;
  }

  point operator*(fuzzy a) const { return point{a * real(), a * imag()}; }
  point operator/(fuzzy a) const { return point{a / real(), a / imag()}; }

  point operator+(point const &a) const {
    return static_cast<std::complex<fuzzy>>(a) +
           static_cast<std::complex<fuzzy>>(*this);
  }
  point operator-(point const &a) const {
    return static_cast<std::complex<fuzzy>>(a) -
           static_cast<std::complex<fuzzy>>(*this);
  }
  point operator*(point const &a) const {
    return static_cast<std::complex<fuzzy>>(a) *
           static_cast<std::complex<fuzzy>>(*this);
  }
  point operator/(point const &a) const {
    return static_cast<std::complex<fuzzy>>(a) *
           static_cast<std::complex<fuzzy>>(*this);
  }
  bool operator==(point const &a) const {
    return static_cast<std::complex<fuzzy>>(a) ==
           static_cast<std::complex<fuzzy>>(*this);
  }
  bool operator!=(point const &a) const {
    return static_cast<std::complex<fuzzy>>(a) !=
           static_cast<std::complex<fuzzy>>(*this);
  }
};

point operator*(fuzzy a, point const &p) {
  return point{a * p.real(), a * p.imag()};
}
point operator/(fuzzy a, point const &p) {
  return point{a / p.real(), a / p.imag()};
}

bool operator<(point const &a, point const &b) {
  return a.real() != b.real() ? a.real() < b.real() : a.imag() < b.imag();
}

using points = std::vector<point>;

struct line : std::pair<point, point> {
  line() = default;
  line(point const &a, point const &b) : std::pair<point, point>(a, b) {}

  explicit operator point() const { return point{second - first}; }

  bool collide(point const &b) const { return abs(first.ccw(second, b)) != 1; }
  bool collide(line const &b) const {
    line l{static_cast<point>(*this), static_cast<point>(b)};
    return !l.collide(point{}) || collide(b.second);
  }
  bool collide(segment const &b) const;

  point proj(point const &p) const {
    point x = {second - first};
    return first + x * (x.dot(p - first) / norm(x));
  }

  point reflection(point const &p) const { return 2.0 * proj(p) - p; }
  point cross_point(line const &b) const {
    fuzzy d1 = static_cast<point>(b).cross({b.first - first});
    fuzzy d2 = static_cast<point>(b).cross(second - first);

    if (d1 == 0 && d2 == 0)
      return first; // same line
    if (d2 == 0)
      throw std::invalid_argument{"No cross point"};
    return first + static_cast<fuzzy::value_type>(d1 / d2) * (second - first);
  }

  fuzzy dist(point const &p) const { return abs(proj(p) - p); }
  fuzzy dist(line const &l) const {
    return collide(l) ? fuzzy{0} : dist(l.first);
  }
  fuzzy dist(segment const &s) const;
  fuzzy dist(circle const &c) const;

  points cross_points(circle const &c) const;
};

using lines = std::vector<line>;

struct segment : line {
  segment(point const &a, point const &b) : line(a, b) {}

  bool collide(segment const &b) const {
    return first.ccw(second, b.first) * first.ccw(second, b.second) <= 0 &&
           b.first.ccw(b.second, first) * b.first.ccw(b.second, second) <= 0;
  }
  bool collide(point const &b) const { return !first.ccw(second, b); }

  fuzzy dist(point const &p) const {
    point r = proj(p);
    if (collide(r)) {
      return abs(r - p);
    }
    return std::min(abs(first - p), abs(second - p));
  }
  fuzzy dist(segment const &s) const {
    if (collide(s)) {
      return 0.0;
    }
    return std::min(std::min(dist(s.first), dist(s.second)),
                    std::min(s.dist(first), s.dist(second)));
  }
  fuzzy dist(circle const &c) const;
};

bool line::collide(segment const &b) const {
  point p{b.first - first};
  return p.cross(b.first - first) * p.cross(b.second - first) < 0;
}

fuzzy line::dist(segment const &s) const {
  return collide(s) ? fuzzy{0} : std::min(dist(s.first), dist(s.second));
}

struct circle {
  point center;
  fuzzy radius;

  fuzzy dist(circle const &c) const {
    fuzzy d = abs(center - c.center);
    return d < fuzzy{abs(radius - c.radius)}
               ? std::max(static_cast<fuzzy::value_type>(d - radius - c.radius),
                          0.0)
               : abs(radius - c.radius) - d;
  }

  points cross_points(circle const &c) const {
    points result{2};
    point ab = c.center - center;
    fuzzy d = abs(ab);
    fuzzy crL = (norm(ab) + radius * radius - c.radius * c.radius) / (2 * d);
    if (d == 0 || radius < fuzzy{abs(crL)})
      return result;

    point abN = ab * point{0, sqrt(radius * radius - crL * crL) / d};
    point cp = center + crL / d * ab;
    result.emplace_back(cp + abN);
    if (norm(abN) != 0)
      result.emplace_back(cp - abN);
    return result;
  }

  points tangent_points(point p) {
    points result{2};
    fuzzy sine = radius / abs(p - center);
    if (sine <= 1)
      return result;

    fuzzy t = M_PI_2 - asin(std::min(sine, fuzzy{1.0}));
    result.emplace_back(center + (p - center) * point(std::polar(sine, t)));
    if (sine != 1)
      result.emplace_back(center + (p - center) * point(std::polar(sine, -t)));
    return result;
  }

  lines common_tangent_lines(circle const &c) {
    lines result;
    fuzzy d = abs(c.center - center);
    REPI(i, 2) {
      fuzzy sine = (radius - (1 - i * 2) * c.radius) / d;
      if (sine * sine <= 1)
        break;

      fuzzy cos = std::sqrt(std::max(1 - sine * sine, 0.0));
      REPI(j, 2) {
        point n = point{sine, (1 - j * 2) * cos} * point(c.center - center) / d;
        result.emplace_back(center + radius * n,
                            c.center + (1 - i * 2) * c.radius * n);
        if (cos < 0)
          break;
      }
    }
    return result;
  }
};

fuzzy line::dist(circle const &c) const {
  return std::max(static_cast<fuzzy::value_type>(dist(c.center) - c.radius),
                  0.0);
}

points line::cross_points(circle const &c) const {
  points result{2};
  point ft{proj(c.center)}, p{static_cast<point>(*this)};
  auto r = c.radius;
  if (r * r >= norm(ft - c.center))
    return result;

  point dir =
      std::sqrt(std::max(
          static_cast<fuzzy::value_type>(r * r - norm(ft - c.center)), 0.0)) /
      norm(p) * (p);
  result.emplace_back(ft + dir);
  if (r * r != norm(ft - c.center)) {
    result.emplace_back(ft - dir);
  }
  return result;
}

fuzzy segment::dist(circle const &c) const {
  fuzzy dSqr1 = norm(c.center - first), dSqr2 = norm(c.center - second);
  auto r = c.radius;
  if (dSqr1 < r * r ^ dSqr2 < r * r)
    return 0;
  if (dSqr1 < r * r & dSqr2 < r * r)
    return r - sqrt(std::max(dSqr1, dSqr2));
  return std::max(static_cast<fuzzy::value_type>(dist(c.center) - r), 0.0);
}

// 三角形の外心
// 条件: a,b,cは同一線上でない
point circumcenter(point const &a, point const &b, point const &c) {
  point x = (a - c) * 0.5;
  point y = (b - c) * 0.5;
  return c + line{x, x * point{1, 1}}.cross_point({y, y * point{1, 1}});
}

// 点aと点bを通る半径がrの円において、中心になりえる点
points circles_points_radius(point const &a, point const &b, fuzzy r) {
  points result;
  point abH = (b - a) * 0.5;
  fuzzy d = abs(abH);
  if (d == 0 || d > r)
    return result; // Do round with (d <= r) condition if needed

  fuzzy dN = sqrt(r * r - d * d); // Do max(r*r - d*d, 0) if needed
  point n = abH * point(0, 1) * (dN / d);
  result.emplace_back(a + abH + n);
  if (dN > 0)
    result.emplace_back(a + abH - n);

  return result;
}

// 点aと点bを通り直線lに接する円において、中心になりえる点
points circles_points_tangent(point const &a, point const &b, line const &l) {
  point n = static_cast<point>(l) * point{0, 1};
  point m = (b - a) * point{0, 0.5};
  fuzzy nm = n.dot(m);
  fuzzy rC = point((a + b) * 0.5 - l.first).dot(n);
  fuzzy qa = norm(n) * norm(m) - nm * nm;
  fuzzy qb = -rC * nm;
  fuzzy qc = norm(n) * norm(m) - rC * rC;
  fuzzy qd = qb * qb - qa * qc;

  points result{2};
  if (qd < -0) {
    return result;
  }
  if (qa == 0) {
    if (qb != 0) {
      result.emplace_back((a + b) * 0.5 - m * (qc / qb / 2));
    }
    return result;
  }

  fuzzy t = -qb / qa;
  result.emplace_back((a + b) * 0.5 + m * (t + sqrt(std::max(qd, {0})) / qa));
  if (qd > 0) {
    result.emplace_back((a + b) * 0.5 + m * (t - sqrt(std::max(qd, {0})) / qa));
  }
  return result;
}

// 点集合を含む最小の円の中心
point min_enclosing_circle(points const &ps) {
  point c;
  double move = 0.5;
  REPI(i, 39) { // 2^(-39-1) ≒ 0.9e-12
    REPI(t, 50) {
      fuzzy max = 0;
      int k = 0;
      REPI(j, ps.size()) {
        if (max < norm(ps[j] - c)) {
          max = norm(ps[j] - c);
          k = j;
        }
      }
      c += (ps[k] - c) * move;
    }
    move /= 2;
  }
  return c;
}

//
// Polygon processes
//

// Convex hull
points convex_hull(points const &ps) {
  auto n = ps.size();
  int k = 0;
  points ch{2 * n};
  for (int i = 0; i < n;) { // lower-hull
    while (k >= 2 && ch[k - 2].ccw(ch[k - 1], ps[i]) <=
                         0) { // Do (== -1) if include extra vertex
      --k;
    }
    ch[k] = ps[i];
    ++k;
    ++i;
  }
  for (int i = n - 2, t = k + 1; i >= 0;) { // upper-hull
    while (k >= t && ch[k - 2].ccw(ch[k - 1], ps[i]) <= 0) {
      --k;
    }
    ch[k] = ps[i];
    ++k;
    --i;
  }
  ch.resize(k - 1);
  return ch;
}

points convex_hull(points &ps) {
  sort(ps.begin(), ps.end());
  return convex_hull(const_cast<points const &>(ps));
}

bool is_ccw_convex(points const &ps) {
  auto n = ps.size();
  REPI(i, n) {
    if (ps[i].ccw(ps[(i + 1) % n], ps[(i + 2) % n]) ==
        -1) { // Do (!= 1) if disallow Degeneracy
      return false;
    }
  }
  return true;
}

// is in a convex polygon O(log n)
// It returns
//   0: outer
//   1: inner
//   2: on the line segment
int in_convex(point const &p, points const &ps) {
  auto n = ps.size();
  point g = (ps[0] + ps[n / 3] + ps[n * 2 / 3]) / 3.0;
  if (g == p)
    return 1;

  point gp = p - g;

  size_t l = 0, r = n;
  while (l + 1 < r) { // binary search
    int mid = (l + r) / 2;
    point gl = ps[l] - g;
    point gm = ps[mid] - g;
    fuzzy glgp = gl.cross(gp);
    fuzzy gmgp = gm.cross(gp);
    if (gl.cross(gm) > 0) {
      if (glgp >= 0 && gmgp <= 0)
        r = mid;
      else
        l = mid;
    } else {
      if (glgp <= 0 && gmgp >= 0)
        l = mid;
      else
        r = mid;
    }
  }
  r %= n;
  fuzzy cr = point(ps[l] - p).cross(ps[r] - p);
  return cr == 0 ? 2 : (cr < 0 ? 0 : 1);
}

// is in a general polygon
// It returns
//   0: outer
//   1: inner
//   2: on the line segment
int in_polygon(point const &p, points const &ps) {
  auto n = ps.size();
  bool in = false;
  REPI(i, n) {
    point a = ps[i] - p;
    point b = ps[(i + 1) % n] - p;
    if (a.cross(b) == 0 && a.dot(b) > 0)
      return 2;
    if (a.imag() > b.imag())
      swap(a, b);
    if ((a.imag() * b.imag() < 0 ||
         (a.imag() * b.imag() < EPS && b.imag() > EPS)) &&
        a.cross(b) > 0)
      in = !in;
  }
  return in ? 1 : 0;
}

// convex polygon clipping
points convex_cut(const points &ps, point const &a1, point const &a2) {
  auto n = ps.size();
  points result;
  REPI(i, n) {
    int ccwc = a1.ccw(a2, ps[i]);
    if (ccwc != -1) {
      result.emplace_back(ps[i]);
    }
    int ccwn = a1.ccw(a2, ps[(i + 1) % n]);
    if (ccwc * ccwn == -1) {
      line l{a1, a2};
      result.emplace_back(l.cross_point(line{ps[i], ps[(i + 1) % n]}));
    }
  }
  return result;
}

// the index of the convex polygon's diameter
std::pair<int, int> convex_diameter_indexes(points const &ps) {
  auto n = ps.size();
  int i = std::distance(ps.begin(), min_element(ps.begin(), ps.end()));
  int j = std::distance(ps.begin(), max_element(ps.begin(), ps.end()));
  int maxI, maxJ;
  fuzzy maxD = 0;
  REPI(ignore, 2 * n) {
    if (maxD < norm(ps[i] - ps[j])) {
      maxD = norm(ps[i] - ps[j]);
      maxI = i;
      maxJ = j;
    }
    if (point(ps[i] - ps[(i + 1) % n]).cross(ps[(j + 1) % n] - ps[j]) <= 0)
      j = (j + 1) % n;
    else
      i = (i + 1) % n;
  }
  return {maxI, maxJ};
}

// the signed area of a polygon
fuzzy area(points const &ps) {
  fuzzy a = 0;
  REPI(i, ps.size()) { a += ps[i].cross(ps[(i + 1) % ps.size()]); }
  return a / 2;
}

// the centroid of polygon
point centroid(points const &ps) {
  auto n = ps.size();
  fuzzy aSum = 0;
  point c;
  REPI(i, n) {
    fuzzy a = ps[i].cross(ps[(i + 1) % n]);
    aSum += a;
    c += (ps[i] + ps[(i + 1) % n]) * a;
  }
  return 1 / aSum / 3 * c;
}

// the voronoi cell
points voronoi_cell(point const &p, points const &ps, points const &outer) {
  points cl{outer};
  REPI(i, ps.size()) {
    if (norm(ps[i] - p) == 0)
      continue;
    point h = (p + ps[i]) * 0.5;
    cl = convex_cut(cl, h, h + (ps[i] - h) * point{0, 1});
  }
  return cl;
}

//
// Geometry graph theory
//

struct Edge {
  int from, to;
  fuzzy cost;
  Edge(int from, int to, fuzzy cost) : from(from), to(to), cost(cost) {}
  Edge(Edge const &old) : from(old.from), to(old.to), cost(old.cost) {}
};

struct Graph {
  int _n;
  std::vector<std::vector<Edge>> edges;
  Graph(int n) : _n(n), edges(n) {}

  void addEdge(Edge const &e) {
    if (!(0 <= e.from && e.from < _n && 0 <= e.to && e.to < _n)) {
      throw std::out_of_range{""};
    }
    edges[e.from].emplace_back(e);
    edges[e.to].emplace_back(e.to, e.from, e.cost);
  }
};

// 線分アレンジメント
Graph segment_arrangement(std::vector<segment> const &segs, points const &ps) {
  auto n = segs.size();
  auto m = ps.size();
  Graph gr(m);
  std::vector<std::pair<fuzzy, int>> list;
  REPI(i, n) {
    list.clear();
    REPI(j, m) {
      if (segs[i].collide(ps[j]))
        list.emplace_back(norm(segs[i].first - ps[j]), j);
    }
    sort(list.begin(), list.end());
    REPI(j, list.size() - 1) {
      int a = list[j].second;
      int b = list[j + 1].second;
      gr.addEdge(Edge(a, b, abs(ps[a] - ps[b])));
    }
  }
  return gr;
}

Graph segment_arrangement(std::vector<segment> const &segs, points &ps) {
  auto n = segs.size();
  REPI(i, n) {
    ps.emplace_back(segs[i].first);
    ps.emplace_back(segs[i].second);
    REPI(j, i) {
      if (segs[i].collide(segs[j]))
        ps.emplace_back(segs[i].cross_point(segs[j]));
    }
  }
  sort(ps.begin(), ps.end());
  ps.erase(unique(ps.begin(), ps.end()), ps.end());
  return segment_arrangement(segs, const_cast<points const &>(ps));
}

inline void visibility_graph_impl(const points &ps,
                                  const std::vector<points> &objs, size_t n,
                                  Graph &gr, size_t i, size_t j) {
  point a = ps[i], b = ps[j];
  if (norm(a - b) != 0) {
    FOR(const &obj, objs) {
      int inStA = in_convex(a, obj);
      int inStB = in_convex(b, obj);
      if ((inStA ^ inStB) % 2 ||
          (inStA * inStB != 1 && in_convex((a + b) * 0.5, obj) == 1))
        return;
      REPI(l, obj.size()) {
        point cur = obj[l];
        point next = obj[(l + 1) % obj.size()];
        if (segment{a, b}.collide(segment{cur, next}) &&
            !segment{cur, next}.collide(a) && !segment{cur, next}.collide(b))
          return;
      }
    }
  }
  gr.addEdge(Edge(i, j, abs(a - b)));
}

// 可視グラフ（点集合から見える位置へ辺を張ったグラフ）
Graph visibility_graph(const points &ps, const std::vector<points> &objs) {
  auto n = ps.size();
  Graph gr(n);
  REPI(i, n) {
    REPI(j, i) { visibility_graph_impl(ps, objs, n, gr, i, j); }
  }
  return gr;
}

//
// Others
//

// 重複する線分を併合する
lines mergeSegments(lines &segs) {
  auto n = segs.size();
  REPI(i, n) {
    if (segs[i].second < segs[i].first)
      swap(segs[i].second, segs[i].first);
  }

  REPI(i, n) {
    REPI(j, i) {
      line &l1 = segs[i], &l2 = segs[j];
      if (point{l1.second - l1.first}.cross(l2.second - l2.first) == 0 &&
          l1.collide(l2.first) && l1.first.ccw(l1.second, l2.second) != 2 &&
          l2.first.ccw(l2.second, l1.second) != 2) {
        segs[j] =
            line(std::min(l1.first, l2.first), std::max(l1.second, l2.second));
        segs[i--] = segs[--n];
        break;
      }
    }
  }
  segs.resize(n);
  return segs;
}

#endif // E_SMART_LIB_GEOMETRY

lib_io.hpp

#ifndef E_SMART_LIB_IO
#define E_SMART_LIB_IO

#include "lib_defs.hpp"

template <class T> T speedy_decimal_read() {
  bool sign = false;
  unsigned char c = '\0';
  do {
    c = getchar();
    if (c == '-')
      sign = true;
  } while (!('0' <= c && c <= '9'));
  T num = 0;
  while ('0' <= c && c <= '9') {
    num = num * 10 + (c - '0');
    c = getchar();
  }
  return sign ? -num : num;
}

#endif // E_SMART_LIB_IO

lib_types.hpp

#ifndef E_SMART_LIB_TYPES
#define E_SMART_LIB_TYPES

#include <cmath>

static constexpr long double EPS = 1e-9;

struct fuzzy {
  using value_type = double;
  value_type v = 0;

  fuzzy(value_type const &value) : v(value) {}
  fuzzy(fuzzy const &src) : v(src.v) {}
  operator value_type() const { return v; }

  template <class T> fuzzy &operator=(T const &r) {
    v = r;
    return *this;
  }
  template <class T> fuzzy operator+(T const &r) const { return v + r; }
  template <class T> fuzzy operator-(T const &r) const { return v - r; }
  template <class T> fuzzy operator*(T const &r) const { return v * r; }
  template <class T> fuzzy operator/(T const &r) const { return v / r; }
  fuzzy operator+() const { return +v; }
  fuzzy operator-() const { return -v; }
  template <class T> fuzzy &operator+=(T const &r) {
    v += r;
    return *this;
  }
  template <class T> fuzzy &operator-=(T const &r) {
    v -= r;
    return *this;
  }
  template <class T> fuzzy &operator*=(T const &r) {
    v *= r;
    return *this;
  }
  template <class T> fuzzy &operator/=(T const &r) {
    v /= r;
    return *this;
  }
  fuzzy &operator++() {
    ++v;
    return *this;
  }
  fuzzy &operator--() {
    --v;
    return *this;
  }
  fuzzy operator++(int) {
    auto self = *this;
    ++v;
    return self;
  }
  fuzzy operator--(int) {
    auto self = *this;
    --v;
    return self;
  }
  template <class T> bool operator<(fuzzy const &r) { return v < r + EPS; }
  template <class T> bool operator>(fuzzy const &r) { return v + EPS > r; }
  template <class T> bool operator<=(fuzzy const &r) { return v + EPS >= r; }
  template <class T> bool operator>=(fuzzy const &r) { return v <= r + EPS; }
  template <class T> bool operator==(fuzzy const &r) {
    return diff_abs(v, r) < EPS;
  }
  template <class T> bool operator!=(fuzzy const &r) {
    return diff_abs(v, r) >= EPS;
  }
};

bool signbit(fuzzy const &x) { return signbit(x.v); }

struct modl {
  using value_type = long;
  value_type v;
  value_type mod;

  modl(value_type const &modulo) : v(0), mod(modulo) {}
  modl(value_type const &value, value_type const &modulo)
      : v(value), mod(modulo) {}
  operator value_type() const { return v; }
  modl &operator=(modl const &r) {
    v = r.v;
    return *this;
  }
  template <class T> modl &operator=(T const &r) {
    v = r;
    return *this;
  }

  modl operator+() const { return {+v, mod}; }
  modl operator-() const { return {-v, mod}; }
  template <class T> modl operator+(T const &r) const {
    return {static_cast<long long>(v) + r % mod, mod};
  }
  template <class T> modl operator-(T const &r) const {
    return {static_cast<long long>(v) - r % mod, mod};
  }
  template <class T> modl operator*(T const &r) const {
    return {static_cast<long long>(v) * r % mod, mod};
  }
  template <class T> modl operator/(T const &r) const {
    return {static_cast<long long>(v) / r % mod, mod};
  }
  template <class T> modl operator%(T const &r) const {
    return {static_cast<long long>(v) % r % mod, mod};
  }
  template <class T> modl &operator+=(T const &r) {
    v += r;
    return *this;
  }
  template <class T> modl &operator-=(T const &r) {
    v -= r;
    return *this;
  }
  template <class T> modl &operator*=(T const &r) {
    v *= r;
    return *this;
  }
  template <class T> modl &operator/=(T const &r) {
    v /= r;
    return *this;
  }
  template <class T> modl &operator%=(T const &r) {
    v %= r;
    return *this;
  }
  modl &operator++() {
    ++v;
    return *this;
  }
  modl &operator--() {
    --v;
    return *this;
  }
};

#endif // E_SMART_LIB_TYPES