【模板】【计几】圆的反演

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 #define eps 1e-6
 4 int sgn(double x){
 5     if( fabs(x) < eps) return 0;
 6     if( x < 0) return -1;
 7     return 1;
 8 }
 9 struct Point{
10     double x,y;
11     Point operator - (const Point& b){
12         return (Point){ x - b.x,y-b.y};
13     }
14     Point operator + (const Point& b){
15         return (Point){x+b.x,y+b.y};
16     }
17     Point operator * (const double& b){
18         return (Point){b * x,b * y};
19     }
20     Point Move(double a,double d){
21         return (Point){ x + d * cos(a),y + d * sin(a)};
22     }
23 };
24 struct Circle{
25     Point o;
26     double r;
27 }c[3],c0,ansc[3];
28 int tot;
29 double cross(Point a,Point b,Point c){
30     return (b.x - a.x) * (c.y - a.y) - (c.x - a.x)*(b.y - a.y);
31 }
32 double dot(Point a,Point b,Point c){
33     return (b.x - a.x) * (c.x - a.x) + (b.y - a.y) * (c.y - a.y);
34 }
35 double dis(Point a,Point b){
36     return sqrt( (a.x - b.x) * (a.x - b.x) + (a.y - b.y)*(a.y - b.y));
37 }
38 Point Point_Inver(Circle c0,Point P){
39     Point OP = P - c0.o;
40     double len = dis(c0.o,P);
41     len = len*len;
42     return c0.o + OP*( c0.r * c0.r / len );
43 }
44 Circle Circle_Inver(Circle c0,Circle a){
45     Circle res;
46     Point OA = a.o - c0.o;
47     double len = dis(a.o,c0.o);
48     Point up = c0.o + OA * ( ( len + a.r) / len );
49     Point down = c0.o + OA *( (len - a.r) / len );
50     up = Point_Inver(c0,up);
51     down = Point_Inver(c0,down);
52     res.o = (up+down) * 0.5;
53     res.r = dis(up,down) * 0.5;
54     return res;
55 }
56 Circle Line_Inver(Circle c0,Point a,Point b){
57     Circle res;
58     double d = fabs( cross(a,c0.o,b) / dis(a,b));
59     res.r = c0.r * c0.r / (2.0 * d);
60 
61     double len = dot(a,b,c0.o) / dis(a,b);
62     Point AB = b - a;
63     Point c = a + AB * (len/dis(a,b));
64     Point CO = c - c0.o;
65     res.o = c0.o + CO * (res.r/d);
66 
67     //double len = dis(a,c[1].o);
68     //res.o = c0.o + (a-c[1].o) * (res.r/len);
69     return res;
70 }
71 void solve(){
72     for(int i = 1;i<=2;++i) c[i] = Circle_Inver(c0,c[i]);
73     if( c[1].r < c[2].r - eps) swap(c[1],c[2]);
74     Point v = c[2].o - c[1].o;
75     double a1 = atan2(v.y,v.x);
76     double a2 = acos( (c[1].r - c[2].r) / dis(c[1].o,c[2].o));
77     Point p1 = c[1].o.Move(a1 + a2,c[1].r);
78     Point p2 = c[2].o.Move(a1 + a2,c[2].r);
79     //cerr<<p1.x<<" "<<p1.y<<" "<<p2.x<<" "<<p2.y<<endl;
80     if( sgn(cross(c[1].o,p1,p2)) == sgn(cross(c0.o,p1,p2)) ) ansc[++tot] = Line_Inver(c0,p1,p2);
81     p1 = c[1].o.Move(a1-a2,c[1].r);
82     p2 = c[2].o.Move(a1-a2,c[2].r);
83     //cerr<<p1.x<<" "<<p1.y<<" "<<p2.x<<" "<<p2.y<<endl;
84     if( sgn(cross(c[1].o,p1,p2)) == sgn(cross(c0.o,p1,p2)) ) ansc[++tot] = Line_Inver(c0,p1,p2);
85 }
86 int main(){
87     c0.r = 10.0;
88     int T; scanf("%d",&T);
89     while(T--){
90         tot = 0;
91         for(int i = 1;i<=2;++i) scanf("%lf %lf %lf",&c[i].o.x,&c[i].o.y,&c[i].r);
92         scanf("%lf %lf",&c0.o.x,&c0.o.y);
93         solve();
94         printf("%d\n",tot);
95         for(int i = 1;i<=tot;++i) printf("%.8f %.8f %.8f\n",ansc[i].o.x,ansc[i].o.y,ansc[i].r);
96     }
97     return 0;
98 }