QOJ.ac

QOJ

IDProblemSubmitterResultTimeMemoryLanguageFile sizeSubmit timeJudge time
#293969#7783. Military ManeuverckisekiWA 10ms10772kbC++2013.1kb2023-12-30 00:49:362023-12-30 00:49:37

Judging History

你现在查看的是最新测评结果

  • [2023-12-30 00:49:37]
  • 评测
  • 测评结果:WA
  • 用时:10ms
  • 内存:10772kb
  • [2023-12-30 00:49:36]
  • 提交

answer

// An AC a day keeps the doctor away.
#include <bits/stdc++.h>

#ifdef local
#define safe std::cerr<<__PRETTY_FUNCTION__<<" line "<<__LINE__<<" safe\n"
#define debug(args...) qqbx(#args, args)
#define orange(args...) danb(#args, args)
using std::cerr;
template <typename ...T> void qqbx(const char *s, T ...args) {
    cerr << "\e[1;32m(" << s << ") = (";
    int cnt = sizeof...(T);
    (..., (cerr << args << (--cnt ? ", " : ")\e[0m\n")));
}
template <typename I> void danb(const char *s, I L, I R) {
    cerr << "\e[1;32m[ " << s << " ] = [ ";
    for (int f = 0; L != R; ++L) cerr << (f++ ? ", " : "") << *L;
    cerr << " ]\e[0m\n";
}
#else
#define safe ((void)0)
#define debug(...) ((void)0)
#define orange(...) ((void)0)
#endif // local
#define all(v) begin(v),end(v)

using namespace std;

namespace geometry {
#define IM imag
#define RE real
using lld = int64_t;
using llf = long double;
using i128 = __int128;
using PT = std::complex<lld>;
using PTF = std::complex<llf>;
using P = PT;

PTF toPTF(PT p) { return PTF{(llf)RE(p), (llf)IM(p)}; }

int sgn(int x) { return (x > 0) - (x < 0); }
int sgn(lld x) { return (x > 0) - (x < 0); }
int sgn(i128 x) { return (x > 0) - (x < 0); }
const llf eps = 1e-9;
int sgn(llf x) { return (x > eps) - (x < -eps); }
lld dot(const P& a, const P& b) { return RE(conj(a) * b); }
lld cross(const P& a, const P& b) { return IM(conj(a) * b); }
llf cross(PTF a, PTF b) {
    return IM(conj(a) * b);
}
int ori(const P& a, const P& b, const P& c) {
  return sgn(cross(b - a, c - a));
}

template <typename V> llf area(const V & pt) {
  llf ret = 0;
  for (int i = 1; i + 1 < (int)pt.size(); i++)
    ret += cross(pt[i] - pt[0], pt[i+1] - pt[0]);
  return ret / 2.0;
}

P rot90(P p) { return P{-IM(p), RE(p)}; }

P readpt() {
    int x, y;
    cin >> x >> y;
    return P(x, y);
}

#define x real()
#define y imag()
PTF getCircum(P a, P b, P c){
    /* lld c1 = dot(a + b, a - b); */
    /* lld c2 = dot(a + c, a - c); */
    /* lld D = cross(a - b, a - c); */
    /* return toPTF((a - b) * c2 + (a - c) * c1) * 0.5il / D; */
    llf a1 = a.x-b.x, b1 = a.y-b.y;
    llf c1 = (a.x+b.x)/2 * a1 + (a.y+b.y)/2 * b1;
    llf a2 = a.x-c.x, b2 = a.y-c.y;
    llf c2 = (a.x+c.x)/2 * a2 + (a.y+c.y)/2 * b2;
    llf X = (c1*b2-b1*c2)/(a1*b2-b1*a2);
    llf Y = (a1*c2-c1*a2)/(a1*b2-b1*a2);
    return PTF(X, Y);
}
#undef x
#undef y

bool cmpxy(const P &a, const P &b) {
    return make_pair(RE(a), IM(a)) < make_pair(RE(b), IM(b));
}
// from NaCl, counterclockwise, be careful of n<=2
vector<P> convex_hull(vector<P> v) {
  sort(all(v), cmpxy); // by X then Y
  if (v[0] == v.back()) return {v[0]};
  int t = 0, s = 1; vector<P> h(v.size() + 1);
  for (int _ = 2; _--; s = t--, reverse(all(v)))
    for (P p : v) {
      while (t>s && ori(p, h[t-1], h[t-2]) >= 0) t--;
      h[t++] = p;
    }
  return h.resize(t), h;
}

namespace delaunay_impl {/*{{{*/
const int maxn = 3000;
const int C = 1e9;
/* A triangulation such that no points will strictly
inside circumcircle of any triangle.
find(root, p) : return a triangle contain given point
add_point : add a point into triangulation
Region of triangle u: iterate each u.e[i].tri,
each points are u.p[(i+1)%3], u.p[(i+2)%3]
Voronoi diagram: for each triangle in `res`,
the bisector of all its edges will split the region. */
#define L(i) ((i)==0 ? 2 : (i)-1)
#define R(i) ((i)==2 ? 0 : (i)+1)
#define FOR for (int i = 0; i < 3; i++)
bool in_cc(const array<P,3> &p, P q) {
    i128 det = 0;
    FOR det += i128(norm(p[i]) - norm(q)) *
      cross(p[R(i)] - q, p[L(i)] - q);
    return det > 0;
}
struct Tri;
struct E {
  Tri *t; int side; E() : t(0), side(0) {}
  E(Tri *t_, int side_) : t(t_), side(side_){}
};
struct Tri {
  array<P,3> p; array<Tri*,3> ch; array<E,3> e;
  Tri(){} Tri(P a, P b, P c) : p{a, b, c}, ch{} {}
  bool has_chd() const { return ch[0] != nullptr; }
  bool contains(P q) const {
    FOR if (ori(p[i], p[R(i)], q) < 0) return false;
    return true;
  }
} pool[maxn * 10], *it;
void link(E a, E b) {
  if (a.t) a.t->e[a.side] = b;
  if (b.t) b.t->e[b.side] = a;
}
struct Trigs {
  Tri *root;
  Trigs() { // should at least contain all points
    root =  // C is recommended to be about 100*MAXC^2
      new(it++) Tri(P(-C, -C), P(C*2, -C), P(-C, C*2));
  }
  void add_point(P p) { add_point(find(p, root), p); }
  static Tri* find(P p, Tri *r) {
    while (r->has_chd()) for (Tri *c: r->ch)
        if (c && c->contains(p)) { r = c; break; }
    return r;
  }
  void add_point(Tri *r, P p) {
    array<Tri*, 3> t; /* split into 3 triangles */
    FOR t[i] = new(it++) Tri(r->p[i], r->p[R(i)], p);
    FOR link(E(t[i], 0), E(t[R(i)], 1));
    FOR link(E(t[i], 2), r->e[L(i)]);
    r->ch = t;
    FOR flip(t[i], 2);
  }
  void flip(Tri* A, int a) {
    auto [B, b] = A->e[a]; /* flip edge between A,B */
    if (!B || !in_cc(A->p, B->p[b])) return;
    Tri *X = new(it++)Tri(A->p[R(a)],B->p[b],A->p[a]);
    Tri *Y = new(it++)Tri(B->p[R(b)],A->p[a],B->p[b]);
    link(E(X,0), E(Y,0));
    link(E(X,1), A->e[L(a)]); link(E(X,2), B->e[R(b)]);
    link(E(Y,1), B->e[L(b)]); link(E(Y,2), A->e[R(a)]);
    A->ch = B->ch = {X, Y, nullptr};
    flip(X, 1); flip(X, 2); flip(Y, 1); flip(Y, 2);
  }
};
vector<Tri*> res; set<Tri*> vis;
void go(Tri *now) { // store all tri into res
  if (!vis.insert(now).second) return;
  if (!now->has_chd()) return res.push_back(now);
  for (Tri *c: now->ch) if (c) go(c);
}
void build(vector<P> &ps) {
  it = pool; res.clear(); vis.clear();
  shuffle(ps.begin(), ps.end(), mt19937(114514));
  Trigs tr; for (P p: ps) tr.add_point(p);
  go(tr.root); // use `res` afterwards
}
#undef L
#undef R
#undef FOR
}
vector<vector<int>> nearest_delaunay(vector<P> ps) {
    map<pair<lld, lld>, int> mp;
    for (int i = 0; i < (int)ps.size(); i++) {
        auto pp = ps[i];
        mp[{RE(pp), IM(pp)}] = i + 1;
    }
    delaunay_impl::build(ps);
    vector<vector<int>> g(ps.size());
    for (auto t: delaunay_impl::res) {
        for (int j = 0; j < 3; j++) {
            P a = t->p[j], b = t->p[(j+1)%3];
            int ai = mp[{RE(a), IM(a)}];
            int bi = mp[{RE(b), IM(b)}];
            if (!ai || !bi) continue;
            --ai, --bi;
            g[ai].emplace_back(bi);
            g[bi].emplace_back(ai);
        }
    }
    return g;
}/*}}}*/

int quad(P p) {
  return (IM(p) == 0) // use sgn for PTF
    ? (RE(p) < 0 ? 3 : 1) : (IM(p) < 0 ? 0 : 2);
}
int argCmp(P a, P b) {
  // returns 0/+-1, starts from theta = -PI
  int qa = quad(a), qb = quad(b);
  if (qa != qb) return sgn(qa - qb);
  return sgn(cross(b, a));
}

struct Line {
  P st, ed, dir;
  Line (P s, P e) : st(s), ed(e), dir(e - s) {}
  Line() = default;
}; using L = const Line &;
PTF intersect(L A, L B) {
  llf t = cross(B.st - A.st, B.dir) /
    llf(cross(A.dir, B.dir));
  return toPTF(A.st) + toPTF(A.dir) * t; // C^3 / C^2
}

tuple<lld,lld,lld> fraction_intersect(L A, L B) {
  lld u = cross(B.st - A.st, B.dir);
  lld v = cross(A.dir, B.dir);
  P I = A.st * v + A.dir * u;
  return {RE(I), IM(I), v};
}

bool cov(L l, L A, L B) {
  i128 u = cross(B.st-A.st, B.dir);
  i128 v = cross(A.dir, B.dir);
  // ori(l.st, l.ed, A.st + A.dir*(u/v) - l.st) <= 0?
  i128 x = RE(A.dir) * u + RE(A.st - l.st) * v;
  i128 y = IM(A.dir) * u + IM(A.st - l.st) * v;
  return sgn(x*IM(l.dir) - y*RE(l.dir)) * sgn(v) >= 0;
} // x, y are C^3, also sgn<i128> is needed
bool operator<(L a, L b) {
  if (int c = argCmp(a.dir, b.dir)) return c == -1;
  return ori(a.st, a.ed, b.st) < 0;
}
// cross(pt-line.st, line.dir)<=0 <-> pt in half plane
// the half plane is the LHS when going from st to ed
vector<PTF> HPI(vector<Line> &q, bool skip_sort = false) {
  if (!skip_sort) sort(q.begin(), q.end());
  else assert (is_sorted(q.begin(), q.end()));
  int n = static_cast<int>(q.size()), l = 0, r = -1;
  for (int i = 0; i < n; i++) {
    if (i && !argCmp(q[i].dir, q[i-1].dir)) continue;
    while (l < r && cov(q[i], q[r-1], q[r])) --r;
    while (l < r && cov(q[i], q[l], q[l+1])) ++l;
    q[++r] = q[i];
  }
  while (l < r && cov(q[l], q[r-1], q[r])) --r;
  while (l < r && cov(q[r], q[l], q[l+1])) ++l;
  n = r - l + 1; // q[l .. r] are the lines
  if (n <= 1 || !argCmp(q[l].dir, q[r].dir)) {
      q.clear();
      return {};
  }
  vector<PTF> pt(n);
  for (int i = 0; i < n; i++)
    pt[i] = intersect(q[i+l], q[(i+1)%n+l]);
  q = vector(q.begin() + l, q.begin() + r + 1);
  return pt;
} // test @ 2020 Nordic NCPC : BigBrother



}

template <typename T>
bool chmin(T &x, T v) {
    if (x > v)
        return x = v, true;
    return false;
}

/* string to_string(__int128 x) { */
/*     if (x == 0) return "0"; */
/*     if (x < 0) return "-" + to_string(-x); */
/*     string s; */
/*     while (x) { */
/*         s += char('0' + x % 10); */
/*         x /= 10; */
/*     } */
/*     reverse(all(s)); */
/*     return s; */
/* } */


signed main() {
  ios_base::sync_with_stdio(false);
  cin.tie(nullptr);
  constexpr int64_t inf = 1e11;
  constexpr long double eps = 1e-9;

  using geometry::P;
  using geometry::PTF;
  using geometry::lld;
  using geometry::llf;
  using geometry::toPTF;
  using geometry::Line;
  using geometry::rot90;
  using geometry::sgn;

  const vector<Line> frame = [inf] {
    int xl, xr, yl, yr;
    cin >> xl >> yl >> xr >> yr;
    xl *= 2;
    yl *= 2;
    xr *= 2;
    yr *= 2;
    vector<Line> ls;
    ls.emplace_back(P(xr, 0), P(xr, 0) + P(0, 1));
    ls.emplace_back(P(0, yr), P(0, yr) + P(-1, 0));
    ls.emplace_back(P(xl, 0), P(xl, 0) + P(0, -1));
    ls.emplace_back(P(0, yl), P(0, yl) + P(1, 0));
    auto p = HPI(ls);
    assert (p.size() == 4 && ls.size() == 4);
    return ls;
  } ();

  int n;
  cin >> n;
  vector<P> p(n);
  for (int i = 0; i < n; i++) {
    p[i] = geometry::readpt();
  }
  ranges::sort(p, geometry::cmpxy);
  p.erase(unique(all(p)), end(p));
  n = int(p.size());

  const auto h = geometry::convex_hull(p);
  const int m = static_cast<int>(h.size());
  // assert (m >= 2);
  orange(all(p));
  orange(all(h));

  vector<vector<Line>> farthest_bisectors(m);
  vector<vector<Line>> nearest_bisectors(n);

  // m * (m + 4) (HPI without sort)
  for (int i = 0; i < m; i++) {
    vector<Line> ls;
    for (int j = 1; j < m; j++) {
      P mid = h[i] + h[(i + j) % m];
      Line bis(mid, mid + rot90(h[i] - h[(i + j) % m]));
      ls.emplace_back(bis);
    }
    rotate(ls.begin(), min_element(all(ls)), ls.end());
    assert (is_sorted(all(ls)));
    auto &nls = farthest_bisectors[i];
    nls.resize(ls.size() + frame.size());
    merge(all(ls), all(frame), nls.begin());
    HPI(nls, true);
  }

  // n * 16 * log(n) (HPI with sort)
  int tot = 0;
  auto near_adj = geometry::nearest_delaunay(p);
  for (int i = 0; i < n; i++) {
    auto &ls = nearest_bisectors[i];
    ls = frame;
    for (int j: near_adj[i]) {
      P mid = p[i] + p[j];
      Line bis(mid, mid + rot90(p[j] - p[i]));
      ls.emplace_back(bis);
    }
    tot += (int)ls.size();
    HPI(ls);
  }
  assert (tot <= n * 16);

  llf ans = 0;

  vector<Line> ls;
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < m; j++) {
      ls.resize(farthest_bisectors[j].size() + nearest_bisectors[i].size());
      merge(all(farthest_bisectors[j]), all(nearest_bisectors[i]), ls.begin());

      const auto pt = HPI(ls, true);
      orange(all(pt));
      debug(geometry::area(pt));
      debug(i, j);

      const int k = (int)pt.size();
      for (int z = 0; z < 2; z++) {
        for (int t = 0; t < k; t++) {
          PTF O = toPTF(z == 0 ? p[i] : h[j]);
          auto A = pt[t] / llf(2) - O;
          auto B = pt[(t + 1) % k] / llf(2) - O;

          // \int _ l ^ r (ax+b) x^2 dx

          llf cur = 0;

          do {
            const llf l = real(A);
            const llf r = real(B);
            if (abs(l - r) < eps) continue;
            llf a = (imag(B) - imag(A)) / (r - l);
            llf b = imag(A) - l * a;
            cur -= a * (r*r*r*r - l*l*l*l) / 4 + b * (r*r*r - l*l*l) / 3;

            // b * (r - l) = imag(A) * (r - l) - l * (imag(B) - imag(A))
            // auto Z = imag(B) - imag(A);
            // cur -= (r*r+l*l) * (r + l) * Z / 4 + (r*r + l*r + l*l) * (imag(A) * (r - l) - l * Z) / 3;
          } while (false);
          do {
            const llf l = imag(A);
            const llf r = imag(B);
            if (abs(l - r) < eps) continue;
            llf a = (real(B) - real(A)) / (r - l);
            llf b = real(A) - l * a;
            cur += a * (r*r*r*r - l*l*l*l) / 4 + b * (r*r*r - l*l*l) / 3;
            // auto Z = real(B) - real(A);
            // cur += (r*r+l*l) * (r + l) * Z / 4 + (r*r + l*r + l*l) * (real(A) * (r - l) - l * Z) / 3;
          } while (false);


          if (z == 0)
            ans -= cur;
          else
            ans += cur;
        }
      }
    }
  }

  auto f = frame;
  auto A = geometry::area(HPI(f, false)) / 4;
  debug(ans, A, ans / A);
  cout << fixed << setprecision(20);
  ans = max<llf>(0, ans);

  cout << numbers::pi * ans / A << '\n';

}

Details

Tip: Click on the bar to expand more detailed information

Test #1:

score: 100
Accepted
time: 1ms
memory: 7636kb

input:

0 0 2 2
2
3 1
1 3

output:

8.37758040957278164122

result:

ok found '8.3775804', expected '8.3775804', error '0.0000000'

Test #2:

score: 0
Accepted
time: 0ms
memory: 7492kb

input:

0 0 2 2
2
5 1
1 3

output:

37.69911184307751738851

result:

ok found '37.6991118', expected '37.6991118', error '0.0000000'

Test #3:

score: -100
Wrong Answer
time: 10ms
memory: 10772kb

input:

-2911 2151 336 5941
2000
-83 79
-94 47
48 -29
-47 64
84 75
-44 -86
-58 -11
-31 58
20 53
80 -19
-82 74
-60 -26
8 -68
-42 -61
-14 12
-58 -18
92 10
35 -26
71 64
76 89
-80 6
70 4
-96 -99
95 -80
-3 -22
71 -89
-75 17
-35 -82
-59 95
60 48
-74 50
-82 90
-26 5
-75 -31
-45 85
85 14
-70 -57
59 46
55 13
-23 60
...

output:

6630790927.82185489637777209282

result:

wrong answer 1st numbers differ - expected: '6657168.1428534', found: '6630790927.8218546', error = '995.0377724'