fork download
copy 
  1. #include <algorithm>
  2. #include <cmath>
  3. #include <iomanip>
  4. #include <iostream>
  5. #include <vector>
  6. using namespace std;
  7. const int maxn = 5e4 + 5;
  8. const long double pi = acos((long double)(-1));
  9.  
  10. int n, m, beg[maxn], fen[2 * maxn], pos[maxn][2], prv[2 * maxn], nxt[2 * maxn];
  11. long double x, y, a[maxn], b[maxn];
  12. vector<pair<long double, int> > circx, span;
  13.  
  14. void computeCircleIntersections(long double r) {
  15.     circx.clear();
  16.     for (int i = 0; i < n; i++) {
  17.         long double ca = a[i], cb = b[i];
  18.         long double discrim = ca * ca * r * r - cb * cb + r * r;
  19.         if (discrim > 0) {
  20.             long double sq = sqrt(discrim);
  21.             long double x1 = (-sq - ca * cb) / (ca * ca + 1);
  22.             long double x2 = ( sq - ca * cb) / (ca * ca + 1);
  23.             long double y1 = ca * x1 + cb, y2 = ca * x2 + cb;
  24.             long double a1 = atan2(y1, x1), a2 = atan2(y2, x2);
  25.             circx.push_back(make_pair(a1, i));
  26.             circx.push_back(make_pair(a2, i));
  27.         }
  28.     }
  29.     sort(circx.begin(), circx.end());
  30. }
  31.  
  32. void add(int idx, int d) {
  33.     ++idx;
  34.     while (idx < 2 * maxn) {
  35.         fen[idx] += d;
  36.         idx += (idx & -idx);
  37.     }
  38. }
  39.  
  40. int sum(int idx) {
  41.     ++idx;
  42.     int ret = 0;
  43.     while (idx > 0) {
  44.         ret += fen[idx];
  45.         idx -= (idx & -idx);
  46.     }
  47.     return ret;
  48. }
  49.  
  50. long long countCircleIntersections() {
  51.     long long ret = 0;
  52.     for (int i = 0, tot = 0; i < circx.size(); i++) {
  53.         int idx = circx[i].second;
  54.         if (beg[idx] == -1) {
  55.             ++tot;
  56.             beg[idx] = i;
  57.             add(i, 1);
  58.         } else {
  59.             --tot;
  60.             add(beg[idx], -1);
  61.             ret += tot - sum(beg[idx]);
  62.             beg[idx] = -1;
  63.         }
  64.     }
  65.     return ret;
  66. }
  67.  
  68. void initCyclicList() {
  69.     for (int i = 0; i < circx.size(); i++) {
  70.         prv[i] = (i + circx.size() - 1) % circx.size();
  71.         nxt[i] = (i + 1) % circx.size();
  72.     }
  73. }
  74.  
  75. void deleteCyclicListElement(int idx) {
  76.     int pidx = prv[idx];
  77.     int nidx = nxt[idx];
  78.     prv[nidx] = pidx;
  79.     nxt[pidx] = nidx;
  80. }
  81.  
  82. long double intersectionDistance(int idx1, int idx2) {
  83.     --m;
  84.     long double ix = (b[idx1] - b[idx2]) / (a[idx2] - a[idx1]);
  85.     long double iy = a[idx1] * ix + b[idx1];
  86.     long double dist = sqrt(ix * ix + iy * iy);
  87.     return dist;
  88. }
  89.  
  90. long double sumCyclicList(int start, int end, int *it) {
  91.     long double sum = 0;
  92.     for (int idx = it[start]; idx != end; idx = it[idx]) {
  93.         sum += intersectionDistance(circx[start].second, circx[idx].second);
  94.     }
  95.     return sum;
  96. }
  97.  
  98. long double sumCircleIntersections(long double r) {
  99.     computeCircleIntersections(r);
  100.     if (countCircleIntersections() > m) {
  101.         return 0; //degenerate case: a lot of intersections on query point
  102.     }
  103.     for (int i = 0; i < circx.size(); i++) {
  104.         int idx = circx[i].second;
  105.         if (beg[idx] == -1) {
  106.             beg[idx] = i;
  107.             pos[idx][0] = i;
  108.         } else {
  109.             pos[idx][1] = i;
  110.             long double sp = circx[i].first - circx[beg[idx]].first;
  111.             if (2 * pi - sp < sp) {
  112.                 sp = 2 * pi - sp;
  113.                 int tmp = pos[idx][0];
  114.                 pos[idx][0] = pos[idx][1];
  115.                 pos[idx][1] = tmp;
  116.             }
  117.             span.push_back(make_pair(sp, idx));
  118.         }
  119.     }
  120.     sort(span.begin(), span.end());
  121.     initCyclicList();
  122.     long double sum = 0;
  123.     for (int i = 0; i < span.size(); i++) {
  124.         int idx = span[i].second;
  125.         if (pos[idx][0] < pos[idx][1]) {
  126.             sum += sumCyclicList(pos[idx][0], pos[idx][1], nxt);
  127.         } else {
  128.             sum += sumCyclicList(pos[idx][0], pos[idx][1], nxt);
  129.         }
  130.         deleteCyclicListElement(pos[idx][0]);
  131.         deleteCyclicListElement(pos[idx][1]);
  132.     }
  133.     return sum + m * r;
  134. }
  135.  
  136. int main() {
  137.     ios_base::sync_with_stdio(0);
  138.     cin.tie(NULL);
  139.     cin >> n >> x >> y >> m;
  140.     x /= 1000, y /= 1000;
  141.     for (int i = 0; i < n; i++) {
  142.         cin >> a[i] >> b[i];
  143.         a[i] /= 1000, b[i] /= 1000;
  144.         b[i] += a[i] * x - y;
  145.     }
  146.     fill_n(beg, maxn, -1);
  147.     long double low = 0, high = 1e10;
  148.     for (int i = 0; i < 70; i++) {
  149.         long double mid = (low + high) / 2;
  150.         computeCircleIntersections(mid);
  151.         if (countCircleIntersections() < m) {
  152.             low = mid;
  153.         } else {
  154.             high = mid;
  155.         }
  156.     }
  157.     cout << fixed << setprecision(9) << sumCircleIntersections(low) << '\n';
  158.     return 0;
  159. }
  160.  
Success #stdin #stdout 0s 6388KB
comments (?)
stdin
copy 
4
1000 1000 3
1000 0
-1000 0
0 5000
0 -5000
stdout
copy 
14.282170363
https://ideone.com/Di8gnU
language:
C++14 (gcc 8.3)
created:
9 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!