@description@

给定 n 条直线,构造可重点集合 I 为 n 条直线两两的交点集合。

再给定一点 (p, q)。求该点到 I 中的点的距离中前 m 小的距离之和。

原题链接。

@solution@

求前 m 小,而 m 又这么大,不难想到二分答案。

二分答案 r,可以得到一个以 (p, q) 为圆心、以 r 为半径的圆。

我们想要求圆内直线的交点个数 cnt 并判断是否有 cnt >= m。

注意到此时直线在圆中呈现弦的形式。而圆中弦的相交,可以通过将端点转成极角序,然后判断两个区间严格相交(没有包含的情况)来实现。这个直接树状数组随便做。

然而上面只能解决严格在圆内的交点,恰好在圆上的交点,需要判断一个端点连接了多少条不同的直线来求解。

最后我们得到第 m 大的距离在 [l, l + eps] 之间,且这个区间之间只含这一种距离。

怎么求前 m 大之和呢?我们可以一一找出以 l 为半径的这个圆中的直线交点。先找出严格在圆内的:直接作一个二维偏序即可。

此时在圆上的交点个数可能非常多,所以我们不能取 (l + eps) 为半径,不然就会把这些交点算进严格在圆内的。

圆上的交点个数 = (m - 严格在圆内的交点个数)。

@accepted code@

#include <cmath>
#include <cstdio>
#include <vector>
#include <iomanip>
#include <iostream>
#include <algorithm>
using namespace std; #define double long double const int MAXN = 100000;
const double EPS = 1E-16;
const double PI = acos(-1); int dcmp(double x) {
return (fabs(x) <= EPS ? 0 : (x < 0 ? -1 : 1));
} int n, m; double K[MAXN + 5], B[MAXN + 5]; double t[MAXN + 5]; int cnt;
bool cmp(double a, double b) {return dcmp(a - b) < 0;}
bool cmp2(double a, double b) {return dcmp(a - b) == 0;} int tree[MAXN + 5], le[MAXN + 5], ri[MAXN + 5];
int lowbit(int x) {
return x & -x;
}
int sum(int x) {
int ret = 0;
for(int i=x;i;i-=lowbit(i))
ret += tree[i];
return ret;
}
void update(int x, int k) {
for(int i=x;i<=cnt;i+=lowbit(i))
tree[i] += k;
} int tmp[MAXN + 5];
vector<int>v1[MAXN + 5], v2[MAXN + 5];
bool check(double r) {
cnt = 0;
for(int i=1;i<=n;i++) {
double del = K[i]*K[i]*r*r + r*r - B[i]*B[i];
if( dcmp(del) >= 0 ) {
double frac = K[i]*K[i] + 1; del = sqrt(del);
double x1 = -(del + B[i]*K[i])/frac, y1 = -(K[i]*del - B[i])/frac;
t[++cnt] = atan2(y1, x1);
double x2 = -(-del + B[i]*K[i])/frac, y2 = (K[i]*del + B[i])/frac;
t[++cnt] = atan2(y2, x2);
}
}
sort(t + 1, t + cnt + 1, cmp);
cnt = unique(t + 1, t + cnt + 1, cmp2) - t - 1;
for(int i=1;i<=cnt;i++)
tree[i] = tmp[i] = 0, v1[i].clear(), v2[i].clear();
for(int i=1;i<=n;i++) {
double del = K[i]*K[i]*r*r + r*r - B[i]*B[i];
if( dcmp(del) >= 0 ) {
double frac = K[i]*K[i] + 1; del = sqrt(del);
double x1 = -(del + B[i]*K[i])/frac, y1 = -(K[i]*del - B[i])/frac;
int p = lower_bound(t + 1, t + cnt + 1, atan2(y1, x1), cmp) - t;
double x2 = -(-del + B[i]*K[i])/frac, y2 = (K[i]*del + B[i])/frac;
int q = lower_bound(t + 1, t + cnt + 1, atan2(y2, x2), cmp) - t;
le[i] = min(p, q), ri[i] = max(p, q);
v1[le[i]].push_back(i), v2[ri[i]].push_back(i);
tmp[le[i]]++;
if( le[i] != ri[i] ) tmp[ri[i]]++;
}
}
int tot = 0;
for(int i=1;i<=cnt;i++) {
for(int j=0;j<v1[i].size();j++) update(i, 1);
for(int j=0;j<v2[i].size();j++)
update(le[v2[i][j]], -1), tot += sum(i-1) - sum(le[v2[i][j]]);
tot += 1LL*tmp[i]*(tmp[i] - 1)/2;
}
return tot >= m;
} struct node{
int le, ri, id;
friend bool operator < (node a, node b) {return a.ri < b.ri;}
}b[MAXN + 5], c[MAXN + 5];
double ans; int nw; double intersect(int i, int j) {
double x = (B[j] - B[i]) / (K[i] - K[j]);
double y = K[i] * x + B[i];
return sqrt(x*x + y*y);
} void solve(int le, int ri) {
if( le == ri ) return ;
int m = (le + ri) >> 1;
solve(le, m), solve(m + 1, ri); int p = le, q = m + 1, t = le;
while( p <= m && q <= ri ) {
if( b[p].le < b[q].le ) {
for(int i=q;i<=ri;i++) {
if( b[i].le >= b[p].ri ) break;
if( b[p].ri < b[i].ri ) nw++, ans += intersect(b[p].id, b[i].id);
}
c[t++] = b[p++];
}
else c[t++] = b[q++];
}
while( p <= m ) c[t++] = b[p++];
while( q <= ri ) c[t++] = b[q++];
for(int i=le;i<=ri;i++) b[i] = c[i];
} double get(double r) {
cnt = 0;
for(int i=1;i<=n;i++) {
double del = K[i]*K[i]*r*r + r*r - B[i]*B[i];
if( dcmp(del) >= 0 ) {
double frac = K[i]*K[i] + 1; del = sqrt(del);
double x1 = -(del + B[i]*K[i])/frac, y1 = -(K[i]*del - B[i])/frac;
t[++cnt] = atan2(y1, x1);
double x2 = -(-del + B[i]*K[i])/frac, y2 = (K[i]*del + B[i])/frac;
t[++cnt] = atan2(y2, x2);
}
}
sort(t + 1, t + cnt + 1, cmp);
cnt = unique(t + 1, t + cnt + 1, cmp2) - t - 1;
int tot = 0;
for(int i=1;i<=n;i++) {
double del = K[i]*K[i]*r*r + r*r - B[i]*B[i];
if( dcmp(del) >= 0 ) {
double frac = K[i]*K[i] + 1; del = sqrt(del);
double x1 = -(del + B[i]*K[i])/frac, y1 = -(K[i]*del - B[i])/frac;
int p = lower_bound(t + 1, t + cnt + 1, atan2(y1, x1), cmp) - t;
double x2 = -(-del + B[i]*K[i])/frac, y2 = (K[i]*del + B[i])/frac;
int q = lower_bound(t + 1, t + cnt + 1, atan2(y2, x2), cmp) - t;
tot++, b[tot].le = min(p, q), b[tot].ri = max(p, q), b[tot].id = i;
}
}
sort(b + 1, b + tot + 1), solve(1, tot);
return ans + 1.0*(m - nw)*r;
} int read() {
int x = 0, f = 1; char ch = getchar();
while( (ch != '-') && (ch > '9' || ch < '0') ) ch = getchar();
if( ch == '-' ) f = -1, ch = getchar();
while( '0' <= ch && ch <= '9' ) x = 10*x + ch - '0', ch = getchar();
return x * f;
} int main() {
double x, y;
n = read(), x = 1.0*read()/1000, y = 1.0*read()/1000, m = read();
for(int i=1;i<=n;i++) {
K[i] = 1.0*read()/1000, B[i] = 1.0*read()/1000;
B[i] = K[i]*x + B[i] - y;
}
/*
for(int i=1;i<=n;i++) {
for(int j=i+1;j<=n;j++)
printf("%d %d : ", i, j), cout << dist(intersect(l[i], l[j]), pnt) << endl;
cout << i << " " << (pnt - a[i]) * l[i].ab << endl;
}
*/
double le = 0, ri = 2E9 + 5000;
for(int i=0;i<=50;i++) {
double mid = (le + ri) / 2;
if( check(mid) ) ri = mid;
else le = mid;
}
cout << fixed << setprecision(9) << get(le) << endl;
}

@details@

本来思路很简单,结果它先是卡我精度然后又卡我常。。。

写几点做题时遇到的问题:

(1)可以先把直线平移,并把 (p, q) 平移成坐标系原点

(2)不要觉得向量表达就是万能的,这道题恰巧向量表达精度损失非常大。。。反而用直线与圆的方程联立来得更直接。

(3)二分次数不要太多小心 TLE。cin 是真的慢,能换掉就换。

最新文章

  1. ASP.NET Core CORS 简单使用
  2. mysql 5.7 忘记root密码解决方法
  3. 成都app开发:架构一个App需要学会哪些技术呢?
  4. 使用匿名函数在回调函数中正确访问JS循环变量
  5. aptana studio 3汉化方法 及支持jquery的方法
  6. eclipse部署Tomcat6 : The server does not support version 3.0 of the JEE Web module specification
  7. STM32 USB-HID通信移植步骤
  8. C#基础知识梳理索引
  9. Redis开启远程登录连接。
  10. 简单实现服务器/客户端的c代码
  11. React Native学习(二)—— 开始一个项目
  12. MAC MAMP 中安装配置使用 ThinkPHP
  13. IT题库5-并发和并行
  14. Consider defining a bean of type &#39;com.lvjing.dao.DeviceStatusMapper&#39; in your configuration.
  15. Matlab文本处理:提取指定内容
  16. SQL SERVER 2012/ 2014 分页,用 OFFSET,FETCH NEXT改写ROW_NUMBER的用法(转)
  17. C# IEqualityComparer类型参数写法
  18. IIS字体文件添加MIME映射
  19. [转载]strtok函数和strtok_r函数
  20. WebApi接口传参不再困惑(4):传参详解

热门文章

  1. 最小生成树 状压+prim hdu2489
  2. EL表达式用法---查询博客
  3. MySQL/MariaDB随笔一
  4. 【JVM】体系结构及其细节
  5. 前后端分离,如何在前端项目中动态插入后端API基地址?(in docker)
  6. 如何提交本地代码到git仓库
  7. 可以Postman,也可以cURL.进来领略下cURL的独门绝技
  8. Android中如何使用列表对话框
  9. Java实现选择问题
  10. Java实现有理数的循环节