Abstract
ural 1913 Titan Ruins: Old Generators Are Fine Too
implementation 几何 圆交
Source
http://acm.timus.ru/problem.aspx?space=1&num=1913
Solution
主要思路就是判断距离到某两个点是r，另一个点是2r的点。
比赛时候用三分乱搞，大概是被卡了精度还是什么的，没有通过。
赛后用圆交重写了一遍，一开始还是没通过。
然后发现自己的圆交模板是针对面积问题的，单点的情况没有考虑清楚，乱改的时候改错了。大致就是 处理圆覆盖的情况时，在预处理统计时算被覆盖了一次，在算两圆交点的时候又被统计了一次 。还是没弄熟悉圆交的原理所致。
Code
#include <iostream>
#include
<cassert>
#include
<cmath>
#include
<cstdio>
#include
<cstring>
#include
<algorithm>
using
namespace
std;
const
double
EPS = 1e-
8
;
const
double
PI = acos(-
1.0
);
inline
int
sgn(
double
x) {
if
(fabs(x)<EPS)
return
0
;
return
x>
0
?
1
:-
1
;
}
struct
sp {
double
x, y;
double
pa;
int
cnt;
sp() {}
sp(
double
a,
double
b): x(a), y(b) {}
sp(
double
a,
double
b,
double
c,
int
d): x(a), y(b), pa(c), cnt(d) {}
bool
operator
<(
const
sp &rhs)
const
{
if
(sgn(pa-rhs.pa)==
0
)
return
cnt>
rhs.cnt;
return
pa<
rhs.pa;
}
void
read() {scanf(
"
%lf%lf
"
, &x, &
y);}
void
write() {printf(
"
%.10f %.10f\n
"
, x, y);}
}t, s1, s2, p[
5
], e[
20
];
double
r, rad[
5
];
int
cover[
5
];
sp
operator
*(
double
d,
const
sp &
v) {
return
sp(d*v.x, d*
v.y);
}
sp
operator
-(
const
sp &u,
const
sp &
v) {
return
sp(u.x-v.x, u.y-
v.y);
}
sp
operator
+(
const
sp &u,
const
sp &
v) {
return
sp(u.x+v.x, u.y+
v.y);
}
double
dot(
const
sp &u,
const
sp &
v) {
return
u.x*v.x+u.y*
v.y;
}
double
det(
const
sp &u,
const
sp &
v) {
return
u.x*v.y-u.y*
v.x;
}
double
dis(
const
sp &u,
const
sp &
v) {
double
dx = u.x-
v.x;
double
dy = u.y-
v.y;
return
sqrt(dx*dx+dy*
dy);
}
double
dissqr(
const
sp &u,
const
sp &
v) {
double
dx = u.x-
v.x;
double
dy = u.y-
v.y;
return
dx*dx+dy*
dy;
}
void
write(sp u, sp v) {
puts(
"
Now we have enough power
"
);
u.write(); v.write();
}
bool
cirint(
const
sp &u,
double
ru,
const
sp &v,
double
rv, sp &a, sp &
b) {
double
d =
dis(u, v);
if
(sgn(d-(ru+rv))>
0
|| sgn(d-fabs(ru-rv))<=
0
)
return
0
;
sp c
= u-
v;
double
ca, sa, cb, sb, csum, ssum;
ca
= c.x/d, sa = c.y/
d;
cb
= (rv*rv+d*d-ru*ru)/(
2
*rv*d), sb = sqrt(
1
-cb*
cb);
csum
= ca*cb-sa*
sb;
ssum
= sa*cb+sb*
ca;
a
= sp(rv*csum, rv*
ssum);
a
= a+
v;
sb
= -
sb;
csum
= ca*cb-sa*
sb;
ssum
= sa*cb+sb*
ca;
b
= sp(rv*csum, rv*
ssum);
b
= b+
v;
return
1
;
}
bool
cu(
int
N, sp p[],
double
r[], sp &
res) {
int
i, j, k;
memset(cover,
0
,
sizeof
cover);
for
(i =
0
; i < N; ++
i)
for
(j =
0
; j < N; ++
j) {
double
rd = r[i]-
r[j];
if
(i!=j && sgn(rd)>
0
&& sgn(dis(p[i], p[j])-rd)<=
0
)
cover[j]
++
;
}
sp a, b;
for
(i =
0
; i < N; ++
i) {
int
ecnt =
0
;
e[ecnt
++] = sp(p[i].x-r[i], p[i].y, -PI,
1
);
e[ecnt
++] = sp(p[i].x-r[i], p[i].y, PI, -
1
);
for
(j =
0
; j < N; ++
j) {
if
(j==i)
continue
;
if
(cirint(p[i], r[i], p[j], r[j], a, b)) {
e[ecnt
++] = sp(a.x, a.y, atan2(a.y-p[i].y, a.x-p[i].x),
1
);
e[ecnt
++] = sp(b.x, b.y, atan2(b.y-p[i].y, b.x-p[i].x), -
1
);
if
(sgn(e[ecnt-
2
].pa-e[ecnt-
1
].pa)>
0
) {
e[
0
].cnt++
;
e[
1
].cnt--
;
}
}
}
sort(e, e
+
ecnt);
int
cnt = e[
0
].cnt;
for
(j =
1
; j < ecnt; ++
j) {
if
(cover[i]+cnt==
N) {
double
a = (e[j].pa+e[j-
1
].pa)/
2
;
res
= p[i]+sp(r[i]*cos(a), r[i]*
sin(a));
return
1
;
}
cnt
+=
e[j].cnt;
}
}
return
0
;
}
bool
solve(
const
sp &o,
const
sp &u,
const
sp &v,
double
r) {
p[
0
] = o, rad[
0
] = r+
r;
p[
1
] = u; p[
2
] = v; rad[
1
] = rad[
2
] =
r;
sp a, b;
if
(cu(
3
, p, rad, a)) {
b
=
0.5
*(a-
o);
b
= o+
b;
write(a, b);
return
1
;
}
else
return
0
;
}
int
main() {
cin
>>
r;
t.read(); s1.read(); s2.read();
sp a, b, c, d;
if
(cirint(s1, r, s2, r, a, b)) {
if
(solve(t, s1, s2, r))
return
0
;
}
else
{
if
(cirint(s1, r, t, r, a, b) &&
cirint(s2, r, t, r, c, d)) {
write(a, c);
return
0
;
}
else
{
if
(solve(s1, t, s2, r))
return
0
;
if
(solve(s2, t, s1, r))
return
0
;
}
}
puts(
"
Death
"
);
return
0
;
}
顺便贴圆交模板，这回的应该没问题了。
#include <cmath>
#include
<algorithm>
using
namespace
std;
const
double
PI = acos(-
1.0
);
const
double
EPS = 1e-
8
;
int
sgn(
double
d) {
if
(fabs(d)<EPS)
return
0
;
return
d>
0
?
1
:-
1
;
}
struct
sp {
double
x, y;
double
pa;
int
cnt;
sp() {}
sp(
double
a,
double
b): x(a), y(b) {}
sp(
double
a,
double
b,
double
c,
int
d): x(a), y(b), pa(c), cnt(d) {}
bool
operator
<(
const
sp &rhs)
const
{
/*
需要处理单点的情况时
if (sgn(pa-rhs.pa)==0) return cnt>rhs.cnt;
*/
return
pa<
rhs.pa;
}
void
read() {scanf(
"
%lf%lf
"
, &x, &
y);}
void
write() {printf(
"
%lf %lf\n
"
, x, y);}
};
sp
operator
+(
const
sp &u,
const
sp &
v) {
return
sp(u.x+v.x, u.y+
v.y);
}
sp
operator
-(
const
sp &u,
const
sp &
v) {
return
sp(u.x-v.x, u.y-
v.y);
}
double
det(
const
sp &u,
const
sp &
v) {
return
u.x*v.y-v.x*
u.y;
}
double
dir(
const
sp &p,
const
sp &u,
const
sp &
v) {
return
det(u-p, v-
p);
}
double
dis(
const
sp &u,
const
sp &
v) {
double
dx = u.x-
v.x;
double
dy = u.y-
v.y;
return
sqrt(dx*dx+dy*
dy);
}
//
计算两圆交点
//
输入圆心坐标和半径
//
返回两圆是否相交(相切不算)
//
如果相交，交点返回在a和b(从a到b为u的交弧)
bool
cirint(
const
sp &u,
double
ru,
const
sp &v,
double
rv, sp &a, sp &
b) {
double
d =
dis(u, v);
/*
需要处理单点的情况时 覆盖的情况在这里不统计
if (sgn(d-(ru+rv))>0 || sgn(d-fabs(ru-rv))<=0) return 0;
*/
if
(sgn(d-(ru+rv))>=
0
|| sgn(d-fabs(ru-rv))<=
0
)
return
0
;
sp c
= u-
v;
double
ca, sa, cb, sb, csum, ssum;
ca
= c.x/d, sa = c.y/
d;
cb
= (rv*rv+d*d-ru*ru)/(
2
*rv*d), sb = sqrt(
1
-cb*
cb);
csum
= ca*cb-sa*
sb;
ssum
= sa*cb+sb*
ca;
a
= sp(rv*csum, rv*
ssum);
a
= a+
v;
sb
= -
sb;
csum
= ca*cb-sa*
sb;
ssum
= sa*cb+sb*
ca;
b
= sp(rv*csum, rv*
ssum);
b
= b+
v;
return
1
;
}
sp e[
222
];
int
cover[
111
];
//
求圆并
//
输入点数N，圆心数组p，半径数组r，答案数组s
//
s[i]表示被至少i个圆覆盖的面积(最普通的圆并就是s[1])
void
circle_union(
int
N, sp p[],
double
r[],
double
s[]) {
int
i, j, k;
memset(cover,
0
,
sizeof
cover);
for
(i =
0
; i < N; ++
i)
for
(j =
0
; j < N; ++
j) {
double
rd = r[i]-
r[j];
//
覆盖的情况在这里统计
if
(i!=j && sgn(rd)>
0
&& sgn(dis(p[i], p[j])-rd)<=
0
)
cover[j]
++
;
}
for
(i =
0
; i < N; ++
i) {
int
ecnt =
0
;
e[ecnt
++] = sp(p[i].x-r[i], p[i].y, -PI,
1
);
e[ecnt
++] = sp(p[i].x-r[i], p[i].y, PI, -
1
);
for
(j =
0
; j < N; ++
j) {
if
(j==i)
continue
;
if
(cirint(p[i], r[i], p[j], r[j], a, b)) {
e[ecnt
++] = sp(a.x, a.y, atan2(a.y-p[i].y, a.x-p[i].x),
1
);
e[ecnt
++] = sp(b.x, b.y, atan2(b.y-p[i].y, b.x-p[i].x), -
1
);
if
(sgn(e[ecnt-
2
].pa-e[ecnt-
1
].pa)>
0
) {
e[
0
].cnt++
;
e[
1
].cnt--
;
}
}
}
sort(e, e
+
ecnt);
int
cnt = e[
0
].cnt;
for
(j =
1
; j < ecnt; ++
j) {
double
pad = e[j].pa-e[j-
1
].pa;
s[cover[i]
+cnt] += (det(e[j-
1
], e[j])+r[i]*r[i]*(pad-sin(pad)))/
2
;
cnt
+=
e[j].cnt;
}
}
}
