题目链接
题目大意:
给定三个圆,询问是否存在点满足该点与三个圆夹角均相等,若存在多组解返回夹角最大值。
圆外一点到两圆夹角均相等:
即 sina = sinb = r1 / d1 = r2 / d2
即 d1 / d2 = r1 / r2
引出阿波罗尼斯圆定义:
平面上相异两点A B, 且P 与A B在同一平面内, PA / PB = k (k > 0 且 k != 1), 则 P点的轨迹是圆,把他称为阿波罗尼斯圆。
注意到当k等于1时P点的运动轨迹在A,B的中垂线上,那么此题就转化为直线求交、直线与圆求交、圆与圆求交的简单问题了。
给出阿氏圆的推导结论:
设A(x1, y1), B (x2, y2) , k = a.r / b.r
那么有 圆心坐标为
{(k * k * x2 - x1)/ (k * k - 1), (k * k * y2 - y1) / (k * k -1)}
半径向量为
{(x2 - x1), (y2 - y1)} * k / (k * k - 1)
输出前记得判断点是否合法(不在任一圆内即可)
#include <bits/stdc++.h>
using namespace std;
// `计算几何模板`
const double eps = 1e-4;
const double inf = 1e20;
const double pi = acos(-1.0);
const int maxp = 1010;
//`Compares a double to zero`
int sgn(double x){
if(fabs(x) < eps)return 0;
if(x < 0)return -1;
else return 1;
}
//square of a double
inline double sqr(double x){return x*x;}
struct Point{
double x,y;
Point(){}
Point(double _x,double _y){
x = _x;
y = _y;
}
void input(){
scanf("%lf%lf",&x,&y);
}
bool operator == (Point b)const{
return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
}
Point operator -(const Point &b)const{
return Point(x-b.x,y-b.y);
}
//叉积
double operator ^(const Point &b)const{
return x*b.y - y*b.x;
}
//点积
double operator *(const Point &b)const{
return x*b.x + y*b.y;
}
//返回长度
double len(){
return hypot(x,y);//库函数
}
//返回长度的平方
double len2(){
return x*x + y*y;
}
//返回两点的距离
double distance(Point p){
return hypot(x-p.x,y-p.y);
}
Point operator +(const Point &b)const{
return Point(x+b.x,y+b.y);
}
Point operator *(const double &k)const{
return Point(x*k,y*k);
}
Point operator /(const double &k)const{
return Point(x/k,y/k);
}
//`计算pa 和 pb 的夹角`
//`就是求这个点看a,b 所成的夹角`
//`测试 LightOJ1203`
double rad(Point a,Point b){
Point p = *this;
return fabs(atan2( fabs((a-p)^(b-p)),(a-p)*(b-p) ));
}
//`化为长度为r的向量`
Point trunc(double r){
double l = len();
if(!sgn(l))return *this;
r /= l;
return Point(x*r,y*r);
}
//`逆时针旋转90度`
Point rotleft(){
return Point(-y,x);
}
//`顺时针旋转90度`
Point rotright(){
return Point(y,-x);
}
//`绕着p点逆时针旋转angle`
Point rotate(Point p,double angle){
Point v = (*this) - p;
double c = cos(angle), s = sin(angle);
return Point(p.x + v.x*c - v.y*s,p.y + v.x*s + v.y*c);
}
};
struct Line{
Point s,e;
Line(){}
Line(Point _s,Point _e){
s = _s;
e = _e;
}
//`根据一个点和倾斜角angle确定直线,0<=angle<pi`
Line(Point p,double angle){
s = p;
if(sgn(angle-pi/2) == 0){
e = (s + Point(0,1));
}
else{
e = (s + Point(1,tan(angle)));
}
}
//求线段长度
double length(){
return s.distance(e);
}
//`返回直线倾斜角 0<=angle<pi`
double angle(){
double k = atan2(e.y-s.y,e.x-s.x);
if(sgn(k) < 0)k += pi;
if(sgn(k-pi) == 0)k -= pi;
return k;
}
//`点和直线关系`
//`1 在左侧`
//`2 在右侧`
//`3 在直线上`
int relation(Point p){
int c = sgn((p-s)^(e-s));
if(c < 0)return 1;
else if(c > 0)return 2;
else return 3;
}
// 点在线段上的判断
bool pointonseg(Point p){
return sgn((p-s)^(e-s)) == 0 && sgn((p-s)*(p-e)) <= 0;
}
//`两向量平行(对应直线平行或重合)`
bool parallel(Line v){
return sgn((e-s)^(v.e-v.s)) == 0;
}
//`直线和线段相交判断`
//`-*this line -v seg`
//`2 规范相交`
//`1 非规范相交`
//`0 不相交`
int linecrossseg(Line v){
int d1 = sgn((e-s)^(v.s-s));
int d2 = sgn((e-s)^(v.e-s));
if((d1^d2)==-2) return 2;
return (d1==0||d2==0);
}
//`两直线关系`
//`0 平行`
//`1 重合`
//`2 相交`
int linecrossline(Line v){
if((*this).parallel(v))
return v.relation(s)==3;
return 2;
}
//`求两直线的交点`
//`要保证两直线不平行或重合`
Point crosspoint(Line v){
double a1 = (v.e-v.s)^(s-v.s);
double a2 = (v.e-v.s)^(e-v.s);
return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
}
//点到直线的距离
double dispointtoline(Point p){
return fabs((p-s)^(e-s))/length();
}
//点到线段的距离
double dispointtoseg(Point p){
if(s == e) return p.distance(s);
if(sgn((p-s)*(e-s))<0 || sgn((p-e)*(s-e))<0)
return min(p.distance(s),p.distance(e));
return dispointtoline(p);
}
//`返回线段到线段的距离`
//`前提是两线段不相交,相交距离就是0了`
double dissegtoseg(Line v){
return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));
}
//`返回点p在直线上的投影`
Point lineprog(Point p){
return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );
}
};
//圆
struct circle{
Point p;//圆心
double r;//半径
circle(){}
circle(Point _p,double _r){
p = _p;
r = _r;
}
//输入
void input(){
p.input();
scanf("%lf",&r);
}
//`点和圆的关系`
//`0 圆外`
//`1 圆上`
//`2 圆内`
int relation(Point b){
double dst = b.distance(p);
if(sgn(dst-r) < 0)return 2;
else if(sgn(dst-r)==0)return 1;
return 0;
}
//`线段和圆的关系`
//`比较的是圆心到线段的距离和半径的关系`
int relationseg(Line v){
double dst = v.dispointtoseg(p);
if(sgn(dst-r) < 0)return 2;
else if(sgn(dst-r) == 0)return 1;
return 0;
}
//`直线和圆的关系`
//`比较的是圆心到直线的距离和半径的关系`
int relationline(Line v){
double dst = v.dispointtoline(p);
if(sgn(dst-r) < 0)return 2;
else if(sgn(dst-r) == 0)return 1;
return 0;
}
//`两圆的关系`
//`5 相离`
//`4 外切`
//`3 相交`
//`2 内切`
//`1 内含`
//`需要Point的distance`
//`测试:UVA12304`
int relationcircle(circle v){
double d = p.distance(v.p);
if(sgn(d-r-v.r) > 0)return 5;
if(sgn(d-r-v.r) == 0)return 4;
double l = fabs(r-v.r);
if(sgn(d-r-v.r)<0 && sgn(d-l)>0)return 3;
if(sgn(d-l)==0)return 2;
if(sgn(d-l)<0)return 1;
return 0;// int 需要返回值(除条件外)
}
//`求两个圆的交点,返回0表示没有交点,返回1是一个交点,2是两个交点`
//`需要relationcircle`
//`测试:UVA12304`
int pointcrosscircle(circle v,Point &p1,Point &p2){
int rel = relationcircle(v);
if(rel == 1 || rel == 5)return 0;
double d = p.distance(v.p);
double l = (d*d+r*r-v.r*v.r)/(2*d);
double h = sqrt(r*r-l*l);
Point tmp = p + (v.p-p).trunc(l);
p1 = tmp + ((v.p-p).rotleft().trunc(h));
p2 = tmp + ((v.p-p).rotright().trunc(h));
if(rel == 2 || rel == 4)
return 1;
return 2;
}
//`求直线和圆的交点,返回交点个数`
int pointcrossline(Line v,Point &p1,Point &p2){
if(!(*this).relationline(v))return 0;
Point a = v.lineprog(p);
double d = v.dispointtoline(p);
d = sqrt(r*r-d*d);
if(sgn(d) == 0){
p1 = a;
p2 = a;
return 1;
}
p1 = a + (v.e-v.s).trunc(d);
p2 = a - (v.e-v.s).trunc(d);
return 2;
}
};
circle a, b, c;
Line get_l(circle a, circle b){
Line l = {b.p ,a.p};
double ang = l.angle();
if(sgn(ang-pi/2)>=0) ang -= pi/2;
else ang += pi/2;
Point mid = (a.p+b.p)/2;
return Line{mid, ang};
}
circle get_c(circle a, circle b){
double k = a.r/b.r;
circle c;
c.p = (Point){(k*k*b.p.x-a.p.x)/(k*k-1), (k*k*b.p.y-a.p.y)/(k*k-1)};
Point cr = {b.p.x-a.p.x, b.p.y-a.p.y};
cr= cr*k/(k*k-1);
c.r = cr.len();
return c;
}
bool judge(Point k){
if(a.relation(k) || b.relation(k) || c.relation(k)) return false;
return true;
}
void output(Point k){
if(judge(k)) printf("%.5lf %.5lf\n", k.x, k.y);
}
void wk(circle a, circle b, circle c){
auto L = get_l(a, b);
auto A = get_c(a, c), B = get_c(b, c);
Point a1, a2;
int t = A.pointcrosscircle(B, a1, a2);
if(t == 0) return ;
else if(t == 1){
if(L.relation(a1) == 3) output(a1);
}
else{
if(L.relation(a1) == 3 && L.relation(a2) == 3){
if(judge(a1) && judge(a2)){
if(sgn(a1.distance(a.p)-a2.distance(a.p))<0) output(a1);
else output(a2);
}
else output(a1), output(a2);
}
else if(L.relation(a1) == 3){
output(a1);
}
else if(L.relation(a2) == 3){
output(a2);
}
}
}
int main(){
a.input(), b.input(), c.input();
if(sgn(a.r-b.r)==0 && sgn(a.r-c.r)==0){
auto A = get_l(a, b), B = get_l(a, c), C = get_l(b, c);
if(A.linecrossline(B) == 2){
Point p = A.crosspoint(B);
output(p);
}
}
else if(sgn(a.r-b.r)!=0 && sgn(a.r-c.r)!=0 && sgn(b.r-c.r)!=0){
auto A = get_c(a, b), B = get_c(a, c), C = get_c(b, c);
Point a1, a2, c1, c2;
int t = A.pointcrosscircle(B, a1, a2);
if(t == 0) return 0;
else if(t == 1){
int t1 = A.pointcrosscircle(C, c1, c2);
if(t1){
if(a1 == c1 || a1 == c2) output(a1);
}
}
else{
int t1 = A.pointcrosscircle(C, c1, c2);
if(t1 == 0) return 0;
else if(t1 == 1){
if(c1 == a1 || c1 == a2) output(c1);
}
else{
if((a1 == c1 && a2 == c2) || (a1 == c2 && a2 == c1)){
if(judge(a1) && judge(a2)){
if(sgn(a1.distance(a.p)-a2.distance(a.p))<0) output(a1);
else output(a2);
}
else output(a1), output(a2);
}
else if(a1 == c1 || a1 == c2){
output(a1);
}
else if(a2 == c1 || a2 == c2){
output(a2);
}
}
}
}
else{
if(sgn(a.r-b.r)==0){
wk(a, b, c);
}
else if(sgn(a.r-c.r)==0){
wk(a, c, b);
}
else{
wk(b, c, a);
}
}
return 0;
}