hdu 4773 圆的反演

Posted 2020-06-24 very_czy

第一次接触反演算法。

通过反演圆以求得反演后的直线。

圆的相切即为直线和圆的相切，切点是关键。

值得注意的是，不过中心的直线反演后得到不过中心的圆，圆圆反演后得到另一个圆，因为中途可以得到一条直线，所以可以简化计算。

反演半径不可以随意选取，过大会导致精度问题。

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <map>
#include <set>
#include <stack>
#include <queue>
#include <string>
#include <iostream>
#include <algorithm>
using namespace std;

#define MEM(a,b) memset(a,b,sizeof(a))
#define REP(i,n) for(int i=1;i<=(n);++i)
#define REV(i,n) for(int i=(n);i>=1;--i)
#define FOR(i,a,b) for(int i=(a);i<=(b);++i)
#define RFOR(i,a,b) for(int i=(a);i>=(b);--i)
#define getmid(l,r) ((l) + ((r) - (l)) / 2)
#define MP(a,b) make_pair(a,b)

typedef long long ll;
typedef pair<int,int> pii;
const int INF = (1 << 30) - 1;
const double eps = 1e-10;

double add(double a,double b) //有误差的浮点加法
{
if(abs(a + b) < eps * (abs(a) + abs(b))) return 0;
return a + b;
}

struct Point
{
double x,y;
Point(double tx = 0,double ty = 0) : x(tx),y(ty) {}
Point operator + (Point p)
{
return Point(add(x,p.x),add(y,p.y));
}
Point operator - (Point p)
{
return Point(add(x,-p.x),add(y,-p.y));
}
Point operator * (double d)
{
return Point(x * d,y * d);
}
Point operator / (double d)
{
return Point(x / d,y / d);
}
Point Move(double a,double d)//从x正方向开始，逆时针
{
return Point(x + d * cos(a),y + d * sin(a));
}
void Read()
{
scanf("%lf%lf",&x,&y);
}
};

struct Circle
{
Point o;
double r;
Circle(double tx = 0,double ty = 0,double tr = 0) : o(tx,ty),r(tr) {}
void Read()
{
o.Read();
scanf("%lf",&r);
}
void Out()
{
printf("%.8f %.8f %.8f\n",o.x,o.y,r);
}
};

int Sign(double x) //判断x的正负
{
return (x > eps) - (x < -eps);
}

double Cross(Point a,Point b,Point c) //a为顶点的叉积
{
return (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);
}

double Dis(Point a,Point b)//距离
{
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

Circle c[5];
Point P;
int T,tot;
double R;

Circle Inver(Circle c1) //圆反演圆过程
{
Circle res;
double oc1 = Dis(P,c1.o);
double k1 = 1.0 / (oc1 - c1.r);
double k2 = 1.0 / (oc1 + c1.r);
res.r = 0.5 * (k1 - k2) * R * R;
double oc2 = 0.5 * (k1 + k2) * R * R;
res.o = P + (c1.o - P) * oc2 / oc1;
return res;
}

void Mark(Point a,Point b) //记下答案圆，直线反演回圆
{
++tot;
double t = fabs(Cross(a,P,b) / Dis(a,b)); //求出P点到直线的距离
c[tot].r = R * R / (2.0 * t);
double d = Dis(a,c[1].o);
c[tot].o = P + (a - c[1].o) * (c[tot].r / d); //因为向量(a,c[1].o)与公切线垂直，所以可以利用其长度相似出P到c[tot]圆心的距离
}

void Solve()
{
REP(i,2) c[i] = Inver(c[i]); //将已知两圆反演
if(c[1].r < c[2].r) swap(c[1],c[2]); //c[1]圆的半径较大
Point tmp = c[2].o - c[1].o;//基础应该转角多少
double a1 = atan2(tmp.y,tmp.x); //atan2的经典应用，求出直线的倾斜角！
double a2 = acos((c[1].r - c[2].r) / Dis(c[1].o,c[2].o));//两圆心应该倾斜多少
Point P1 = c[1].o.Move(a1 + a2,c[1].r);//计算切点
Point P2 = c[2].o.Move(a1 + a2,c[2].r);//计算切点
if(Sign(Cross(P1,c[1].o,P2)) == Sign(Cross(P1,P,P2))) Mark(P1,P2); //保证P与c[1]圆心在公切线同侧
P1 = c[1].o.Move(a1 - a2,c[1].r);//同样要考虑公切线在下（上）方的情况，计算切点
P2 = c[2].o.Move(a1 - a2,c[2].r);//计算切点
if(Sign(Cross(P1,c[1].o,P2)) == Sign(Cross(P1,P,P2))) Mark(P1,P2); //保证P与c[1]圆心在公切线同侧
}

int main()
{
R = 10.0; //自定义的反演半径，不可过大
scanf("%d",&T);
REP(tt,T)
{
tot = 2;
REP(i,2) c[i].Read();
P.Read();
Solve();
printf("%d\n",tot - 2);
FOR(i,3,tot) c[i].Out();
}
return 0;
}