Luzhiled / comp-geometry

C++ geometry library for ICPC
https://luzhiled.github.io/comp-geometry/
6 stars 1 forks source link

verify 追加 aoj2373 #47

Open Luzhiled opened 2 years ago

Luzhiled commented 2 years ago

提出: https://onlinejudge.u-aizu.ac.jp/services/review.html#ACPC_2022_05_23/6666817

#include <bits/stdc++.h>
using namespace std;

// template {{{
using i32 = int;
using u32 = unsigned int;
using i64 = long long;
using u64 = unsigned long long;

#define  range(i, l, r) for (i64 i = (i64)(l); i < (i64)(r); (i) += 1)
#define rrange(i, l, r) for (i64 i = (i64)(r) - 1; i >= (i64)(l); (i) -= 1)

#define  whole(f, x, ...) ([&](decltype((x)) container) { return (f)(  begin(container),  end(container), ## __VA_ARGS__); })(x)
#define rwhole(f, x, ...) ([&](decltype((x)) container) { return (f)( rbegin(container), rend(container), ## __VA_ARGS__); })(x)

#define debug(x) cerr << "(" << __LINE__ << ")" << #x << ": " << (x) << endl

constexpr i32 inf   = 1001001001;
constexpr i64 infll = 1001001001001001001ll;

constexpr i32 dx[] = {0, -1, 1, 0, -1, 1, -1, 1}; 
constexpr i32 dy[] = {-1, 0, 0, 1, -1, -1, 1, 1};

struct IoSetup { IoSetup(i32 x = 15){ cin.tie(0); ios::sync_with_stdio(0); cout << fixed << setprecision(x); cerr << fixed << setprecision(x); } } iosetup;

template <typename T = i64> T input() { T x; cin >> x; return x; }

template <typename T> ostream &operator<<(ostream &os, vector<T> &v) { range(i, 0, v.size()) { os << v[i] << (i + 1 != (int)v.size() ? " " : ""); } return os; } 
template <typename T> istream &operator>>(istream &is, vector<T> &v) { for (T &in : v) is >> in; return is; }

template <typename T1, typename T2> ostream &operator<<(ostream &os, pair<T1, T2> p) { os << "(" << p.first << ", " << p.second << ")"; return os; }
template <typename T1, typename T2> istream &operator>>(istream &is, pair<T1, T2> &p) { is >> p.first >> p.second; return is; }

template <typename T> vector<T> make_vector(size_t a, T b) { return vector<T>(a, b); }
template <typename... Ts> auto make_vector(size_t a, Ts... ts) { return vector<decltype(make_vector(ts...))>(a, make_vector(ts...)); }

template <typename T1, typename T2> inline bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
template <typename T1, typename T2> inline bool chmin(T1 &a, T2 b) { return a > b && (a = b, true); }
// }}}

#include "src/base.hpp"
#include "src/product.hpp"
using namespace geometry;

real_number calc(const vector< real_number > &rs) {
  // n := rs.size(), 3 \leq n
  // maximize \sum{ rs[i] * rs[(i + 1) % n] * sin(ts[i]) }
  //            s.t. \sum{ ts[i] } - 2PI = 0
  //
  // f(ts[0], ts[1], ..., ts[n-1]) = \sum{ rs[i] * rs[(i+1)%n] * sin(ts[i]) }
  // g(ts[0], ts[1], ..., ts[n-1]) = \sum{ ts[i] } - 2PI
  //
  // Lagrange function L :=
  // L(ts[0], ts[1], ..., ts[n-1]) = f - \lambda g
  //                               = \sum{ rs[i] * rs[(i+1)%n] * sin(ts[i]) } - \lamdba (\sum{ ts[i] } - 2PI )
  //
  // rs[0] * rs[1] * cos(ts[0]) - \lamdba
  //   = ...
  //   = rs[i] * rs[(i+1)%n] * cos(ts[i]) - \lamdba
  //   = ...
  //   = 0 
  //
  // \lamdba = rs[0] * rs[1] * cos(ts[0])
  // => rs[i] * rs[(i+1)%n] * cos(ts[i]) - rs[0] * rs[1] * cos(ts[0]) = 0
  // => rs[i] * rs[(i+1)%n] * cos(ts[i]) = rs[0] * rs[1] * cos(ts[0])
  // => cos(ts[i]) = (rs[0] * rs[1])/(rs[i] * rs[(i+1)%n]) * cos(ts[0])
  //
  // as[i] := (rs[0] * rs[1]) / (rs[i] * rs[(i+1)%n])
  // ts[0] = m, ts[i] = acos(as[i] * cos(ts[0])), and \sum{ ts[i] } = 2PI
  // => binary search
  assert(rs.size() >= 3);
  int n = rs.size();

  auto make_ts = [&](real_number t0) {
    vector< real_number > ts;
    range(i, 0, n) {
      real_number a = (rs[0] * rs[1]) / (rs[i] * rs[(i+1)%n]);
      ts.emplace_back(acosl(clamp(a * cos(t0), real_number(-1), real_number(1))));
    }
    return ts;
  };
  real_number lb = 0, ub = PI;
  range(qi, 0, 200) {
    real_number m = (lb + ub) / 2;
    vector< real_number > ts = make_ts(m);
    real_number theta_sum = whole(accumulate, ts, real_number(0));

    if (theta_sum < 2 * PI) lb = m;
    else ub = m;
  }

  real_number S = 0;
  vector< real_number > ts = make_ts(lb);
  if (not equals(whole(accumulate, ts, real_number(0)), 2 * PI)) return S;

  range(i, 0, n) {
    S += rs[i] * rs[(i+1)%n] * sinl(ts[i]);
  }

  return S / 2;
}

void solve() {
  int n = input();
  vector< real_number > rs(n);
  cin >> rs;

  whole(sort, rs);

  real_number ans = 0;
  range(bit, 0, 1 << n) {
    vector< real_number > subset_r;
    range(i, 0, n) {
      if (~bit & (1 << i)) continue;
      subset_r.emplace_back(rs[i]);
    }

    if (subset_r.size() < 3) continue;

    do {
      chmax(ans, calc(subset_r));
    } while (next_permutation(subset_r.begin() + 1, subset_r.end()));
  }

  cout << ans << endl;
}

signed main() {
  solve();
}