圆的反演学习笔记

Part 1:相关性质:

一、反演的概念

设在平面内给定一点O和常数k(k不等于零),对于平面内任意一点A,确定A′,使A′为直线OA上一点,并且有向线段OA与OA′满足OA·OA′=k,我们称这种变换是以O为反演中心,以k为反演幂的反演变换,简称反演。称A′为A关于O(r)的互为反演点。

 

二、作已知点的反演点的方法

给出反演极O和反演幂k>0,作点A的反演点A′。

令k=r^2,作出反演基圆⊙O(r),

1)若点A在⊙O(r)外,则过点A作圆的切线(两条),两个切点相连与OA连线交点就是点A′。

2)若点A在⊙O(r)内,则把上述过程逆过来:连结OA,过点A作直线垂直于OA,直线与⊙O(r)的交点处的切线的交点就是点A′。

3)若点A在⊙O(r)上,反演点A′就是点A自身。

4)O没有反演点

 

三、圆的反演变换

圆在不同情形下的反演成像:

1. 当圆不经过反演中心,它的反演图形仍旧是个不过反演中心的圆,并且反演中心为这两个互为反形的圆的位似中心:

2. 当圆与反演圆相交,交点是保持不变的;

3. 当圆在反演圆的外面的时候,反演成像位于圆的内部;反之,当圆位于反演圆的内部,反演成像位于圆的外部。

4. 当圆经过反演中心,它的反演图形是一条直线。反之,任意一条不过反演中心的直线,其反演成像是一个经过反演中心的圆。

5. 相切两圆反向任相切,且切点不变,若切点是反演中心,则其反象是两条平行直线;两圆相切,若反演中心在某圆上,则为反形为相切的直线与圆;

 

 

Part 2 例题:

一、HDU 6097

题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=6097

题意:给定圆内或圆上的两点Q、P,两点到圆心的距离相同,求圆上一点D到Q、P点的距离之和最小。

题解:点在圆内不太好操作,将两点以圆心O为反演中心反演到圆外,做P点关于圆的反演点P',OPD与ODP'相似,相似比是|OP| : r。故转化问题成求DQ'+DP'的最小值,如果是任意两点的话就是椭圆与圆相切,故若Q'P'与圆有交点,则答案为Q'P',否则,因为DQ'=DP',则D为 Q'P'中点 与 O 的连线 与圆的交点。

(个人觉得三分也行,但是没试orz...日后再补)

好的,然后根据这题可以知道原来的点P,反演后的点P',反演中心O,反演圆上的点D',所构成的两个三角形有相似关系。

code:

#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<vector>
#include<queue>
#include<stack>
#define ll long long
#define db double
#define pb push_back
#define mem(a, b) memset(a, b, sizeof (a))
#define scani(a) scanf("%d", &a)
#define scand(a) scanf("%lf", &a)
using namespace std;

const db eps = 1e-8;
const db pi = acos(-1);
int sign(db x){ if(fabs(x) < eps) return 0;  return x > 0 ? 1 : -1; }
int cmp(db k1, db k2){ return sign(k1-k2); }
int inmid(db k1, db k2, db k3){ return sign(k2-k1)*sign(k3-k1) <= 0; }//k1在k2、k3之间
int intersect(db l1,db r1,db l2,db r2){// 区间相交判定,区间1在区间2前
    if(l1>r1) swap(l1,r1); if(l2>r2) swap(l2,r2); return cmp(r1,l2)!=-1 && cmp(r2,l1)!=-1;
}

struct point{
	db x, y;
	point(){}
	point(db k1, db k2){ x = k1, y = k2; }
	point operator + (const point &k1) const { return point(x+k1.x, y+k1.y); }//向量加法、点+向量=点
	point operator - (const point &k1) const { return point(x-k1.x, y-k1.y); }//向量减法、点-点=向量
	point operator * (db k1) const { return (point){x*k1, y*k1}; }//向量数乘
    point operator / (db k1) const { return (point){x/k1, y/k1}; }//向量数除
    int operator == (const point &k1) const { return cmp(x,k1.x)==0 && cmp(y,k1.y)==0; }//比较两个点(向量)是否相同
    point turn(db k1){return (point){x*cos(k1)-y*sin(k1), x*sin(k1)+y*cos(k1)};}//逆时针旋转
    point turn90(){return (point){-y, x};}//逆时针旋转90度
    bool operator < (const point k1) const{
        int a = cmp(x, k1.x);
        if(a == -1) return 1; else if(a == 1) return 0; else return cmp(y,k1.y)==-1;
    }//比较两个点(向量)的大小,x越小则点越小;若x相等,则y越小点越小。可以实现按点的坐标排序
    db len(){ return sqrt(x*x+y*y); }//向量模长
    db len2(){ return x*x+y*y; }//向量模长的平方
    point dir(){ return (*this)/(*this).len(); }
    point unit(){ db k = len(); return point(x/k, y/k); }//单位向量
    db angle() { return atan2(y, x); }//向量的极角
    point getdel(){if (sign(x)==-1||(sign(x)==0&&sign(y)==-1)) return (*this)*(-1); else return (*this);}
    //当横坐标为负时,或横坐标为0纵坐标为负时,将点按原点做对称??? 角度是[-π/2, π/2]
    int getp() const { return sign(y)==1 || (sign(y)==0&&sign(x)==-1); }
	//判断点在12象限,或者在x的负半轴上??? 角度是(0, π]
    void scan(){ scanf("%lf%lf", &x, &y); };
    void print(){ printf("%.11f %.11f\n", x, y); }
};
int inmid(point k1, point k2, point k3){ return inmid(k1.x,k2.x,k3.x) && inmid(k1.y,k2.y,k3.y); }//k1 在 [k2,k3] 内
point midpo(point k1, point k2){ return (k1+k2)/2; }//得到两点中点
db dis2(point k1, point k2){ return (k1.x-k2.x)*(k1.x-k2.x) + (k1.y-k2.y)*(k1.y-k2.y); }
db dis(point k1, point k2){ return sqrt(dis2(k1, k2)); }
db cross(point k1, point k2){ return k1.x*k2.y - k1.y*k2.x; }//叉乘
db dot(point k1, point k2){ return k1.x*k2.x + k1.y*k2.y; }//点乘

struct line{// p[0]->p[1]
    point p[2];
    line(){}
    line(point k1,point k2){p[0]=k1; p[1]=k2;}
    point& operator [] (int k){return p[k];}
    int include(point k){ return sign(cross(p[1]-p[0],k-p[0]))>=0; }//点在直线左侧的判定
    point dir(){return p[1]-p[0];}//方向向量
    line push(){ // 向外 ( 左 ) 平移 eps
        point delta=(p[1]-p[0]).turn90().unit()*eps;
        return {p[0]-delta, p[1]-delta};
    }
};

point proj(point q, point k1, point k2){ point k=k2-k1; return k1+k*(dot(q-k1,k)/k.len2()); }// q 到直线 k1,k2 的投影

db displ(point q, point k1, point k2){// 点到直线的距离
    if(k1 == k2) return dis(q, k1);
    return fabs(cross(k2-k1, q-k1)) / (k2-k1).len();
}

struct circle{
    point o; db r;
    circle(){}
    circle(point _o, db _r){ o = _o, r = _r; }
    bool include(point k){ return cmp(dis(o, k), r) <= 0; }//点在圆内判定
};

circle getinvertcir(circle c1, circle c0){//c1关于c0的反演圆
    circle c2;
    db x0 = c0.o.x, y0 = c0.o.y, r0 = c0.r, x1 = c1.o.x, y1 = c1.o.y, r1 = c1.r;
    db d01 = dis(c0.o, c1.o);
    c2.r = 0.5*((1/(d01-r1))-(1/(d01+r1)))*r0*r0;
    db d02 = r0*r0/(d01+r1)+c2.r;
    //db _d02 = r0*r0/(d01-r1)-c2.r;
    c2.o.x = x0 + d02/d01*(x1-x0);
    c2.o.y = y0 + d02/d01*(y1-y0);
    return c2;
}

circle getinvertcir(line k1, circle c0){//直线k1关于c0的反演圆
    point a = proj(c0.o, k1[0], k1[1]);
    db oa = dis(c0.o, a);
    db ob = c0.r*c0.r/oa;
    point v = a-c0.o;  v = v/v.len();
    point b = c0.o+v*ob;
    circle res;
    res.o = midpo(c0.o, b);
    res.r = ob/2;
    return res;
}

line getinvertline(circle c1, circle c0){//c1关于c0的反演直线
    point v = c1.o - c0.o;
    v = v / v.len();
    db d = c0.r*c0.r / (2*c1.r);
    point k1 = v * d + c0.o;
    v = v.turn90();
    point k2 = k1 + v;
    return line(k1, k2);
}

point getinvertpoint(point p, circle c){
    point v = (p-c.o).dir();
    db len = c.r*c.r/dis(c.o, p);
    return c.o+v*len;
}

int main(){
    int T;
    scanf("%d", &T);
    circle c0;
    c0.o = point(0, 0);
    while(T--){
        point p, q;
        scanf("%lf", &c0.r);
        p.scan();
        q.scan();
        db ck = dis(p, c0.o);
        if(ck < eps){
            printf("%.7f\n", 2*c0.r);
            continue;
        }
        p = getinvertpoint(p, c0);
        q = getinvertpoint(q, c0);
        db d = displ(c0.o, p, q);
        db ans;
        if(cmp(d, c0.r) <= 0){
            ans = dis(p, q);
            ans *= c0.r/dis(c0.o, q);
        }
        else{
            db d1 = dis(c0.o, midpo(p, q)) - c0.r;
            db d2 = dis(p, q)/2.0;
            ans = 2*(sqrt(d1*d1+d2*d2));
            ans *= c0.r/dis(c0.o, q);
        }
        printf("%.7f\n", ans);
    }
    return 0;
}

 

二、HDU 6158

题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=6158

题意:如下图,输入n求前n个图中标识的小圆的面积之和

题解:根据圆的反演的性质,以两个大圆的切点为反演中心,将两个大圆反演成两条垂直x轴的平行直线,然后很容易可以求得各个小圆反演后的圆心和半径,再将小圆的反形反演回来,计算面积。当面积的值小到一定程度的时候就可以忽略不计啦,所以不用担心超时。

code:

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<vector>
#include<queue>
#include<stack>
#define ll long long
#define db double
#define pb push_back
#define mem(a, b) memset(a, b, sizeof (a))
using namespace std;

const db eps = 1e-13;
const db pi = acos(-1);
int sign(db x){ if(fabs(x) < eps) return 0;  return x > 0 ? 1 : -1; }
int cmp(db k1, db k2){ return sign(k1-k2); }
int inmid(db k1, db k2, db k3){ return sign(k2-k1)*sign(k3-k1) <= 0; }//k1在k2、k3之间
int intersect(db l1,db r1,db l2,db r2){// 区间相交判定,区间1在区间2前
    if(l1>r1) swap(l1,r1); if(l2>r2) swap(l2,r2); return cmp(r1,l2)!=-1 && cmp(r2,l1)!=-1;
}

struct point{
	db x, y;
	point(){}
	point(db k1, db k2){ x = k1, y = k2; }
	point operator + (const point &k1) const { return point(x+k1.x, y+k1.y); }//向量加法、点+向量=点
	point operator - (const point &k1) const { return point(x-k1.x, y-k1.y); }//向量减法、点-点=向量
	point operator * (db k1) const { return (point){x*k1, y*k1}; }//向量数乘
    point operator / (db k1) const { return (point){x/k1, y/k1}; }//向量数除
    int operator == (const point &k1) const { return cmp(x,k1.x)==0 && cmp(y,k1.y)==0; }//比较两个点(向量)是否相同
    point turn(db k1){return (point){x*cos(k1)-y*sin(k1), x*sin(k1)+y*cos(k1)};}//逆时针旋转
    point turn90(){return (point){-y, x};}//逆时针旋转90度
    bool operator < (const point k1) const{
        int a = cmp(x, k1.x);
        if(a == -1) return 1; else if(a == 1) return 0; else return cmp(y,k1.y)==-1;
    }//比较两个点(向量)的大小,x越小则点越小;若x相等,则y越小点越小。可以实现按点的坐标排序
    db len(){ return sqrt(x*x+y*y); }//向量模长
    db len2(){ return x*x+y*y; }//向量模长的平方
    point dir(){ return (*this)/(*this).len(); }
    point unit(){ db k = len(); return point(x/k, y/k); }//单位向量
    db angle() { return atan2(y, x); }//向量的极角
    point getdel(){if (sign(x)==-1||(sign(x)==0&&sign(y)==-1)) return (*this)*(-1); else return (*this);}
    //当横坐标为负时,或横坐标为0纵坐标为负时,将点按原点做对称??? 角度是[-π/2, π/2]
    int getp() const { return sign(y)==1 || (sign(y)==0&&sign(x)==-1); }
	//判断点在12象限,或者在x的负半轴上??? 角度是(0, π]
    void scan(){ scanf("%lf%lf", &x, &y); };
    void print(){ printf("%.11f %.11f\n", x, y); }
};
int inmid(point k1, point k2, point k3){ return inmid(k1.x,k2.x,k3.x) && inmid(k1.y,k2.y,k3.y); }//k1 在 [k2,k3] 内
point midpo(point k1, point k2){ return (k1+k2)/2; }//得到两点中点
db dis2(point k1, point k2){ return (k1.x-k2.x)*(k1.x-k2.x) + (k1.y-k2.y)*(k1.y-k2.y); }
db dis(point k1, point k2){ return sqrt(dis2(k1, k2)); }
db cross(point k1, point k2){ return k1.x*k2.y - k1.y*k2.x; }//叉乘
db dot(point k1, point k2){ return k1.x*k2.x + k1.y*k2.y; }//点乘

struct line{// p[0]->p[1]
    point p[2];
    line(){}
    line(point k1,point k2){p[0]=k1; p[1]=k2;}
    point& operator [] (int k){return p[k];}
    int include(point k){ return sign(cross(p[1]-p[0],k-p[0]))>=0; }//点在直线左侧的判定
    point dir(){return p[1]-p[0];}//方向向量
    line push(){ // 向外 ( 左 ) 平移 eps
        point delta=(p[1]-p[0]).turn90().unit()*eps;
        return {p[0]-delta, p[1]-delta};
    }
};

point proj(point q, point k1, point k2){ point k=k2-k1; return k1+k*(dot(q-k1,k)/k.len2()); }// q 到直线 k1,k2 的投影

struct circle{
    point o; db r;
    circle(){}
    circle(point _o, db _r){ o = _o, r = _r; }
    bool include(point k){ return cmp(dis(o, k), r) <= 0; }//点在圆内判定
    db getarea(){ return pi*r*r; }
};

circle getinvertcir(circle c1, circle c0){//c1关于c0的反演圆
    circle c2;
    db x0 = c0.o.x, y0 = c0.o.y, r0 = c0.r, x1 = c1.o.x, y1 = c1.o.y, r1 = c1.r;
    db d01 = dis(c0.o, c1.o);
    c2.r = 0.5*((1/(d01-r1))-(1/(d01+r1)))*r0*r0;
    db d02 = r0*r0/(d01+r1)+c2.r;
    //db _d02 = r0*r0/(d01-r1)-c2.r;
    c2.o.x = x0 + d02/d01*(x1-x0);
    c2.o.y = y0 + d02/d01*(y1-y0);
    return c2;
}

circle getinvertcir(line k1, circle c0){//直线k1关于c0的反演圆
    point a = proj(c0.o, k1[0], k1[1]);
    db oa = dis(c0.o, a);
    db ob = c0.r*c0.r/oa;
    point v = a-c0.o;  v = v/v.len();
    point b = c0.o+v*ob;
    circle res;
    res.o = midpo(c0.o, b);
    res.r = ob/2;
    return res;
}

line getinvertline(circle c1, circle c0){//c1关于c0的反演直线
    point v = c1.o - c0.o;
    v = v / v.len();
    db d = c0.r*c0.r / (2*c1.r);
    point k1 = v * d + c0.o;
    v = v.turn90();
    point k2 = k1 + v;
    return line(k1, k2);
}

point getinvertpoint(point p, circle c){//p关于c的反演点
    point v = (p-c.o).dir();
    db len = c.r*c.r/dis(c.o, p);
    return c.o+v*len;
}

int main(){
    int T;
    scanf("%d", &T);
    circle c0;
    c0.o = point(0, 0); c0.r = 20;
    while(T--){
        db r1, r2;
        scanf("%lf%lf", &r1, &r2);
        int n;
        scanf("%d", &n);
        if(r1 > r2) swap(r1, r2);
        db d1 = c0.r*c0.r/(2*r1);
        db d2 = c0.r*c0.r/(2*r2);
        point p1 = getinvertpoint(point(2*r1, 0), c0);
        point p2 = getinvertpoint(point(2*r2, 0), c0);
        db d = dis(p1, p2);
        circle c1;
        c1.o = midpo(p1, p2);
        c1.r = d/2;
        db sum = 0, area = 0;
        area = getinvertcir(c1, c0).getarea();
        sum += area;
        for(int i = 2; i <= n; i++){
            if(i%2 == 0){
                c1.o.y += d;
            }
            area = getinvertcir(c1, c0).getarea();
            sum += area;
            if(sign(area) == 0) break;
        }
        printf("%.5f\n", sum);
    }
    return 0;
}

 

三、HDU 4773

题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=4773

题意:给定一个点和两个不经过改点且无公共点的两个圆,求所有经过给定点且与两圆外切的圆的圆心坐标和半径。

题解:感觉很明显得可以看出是圆的反演,将两个圆以给定点为反演中心画出反形,求反形的公切线,求出所有切线的反演圆,判断是否与原来的两个圆外切。

code:

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<vector>
#include<queue>
#include<stack>
#define ll long long
#define db double
#define pb push_back
#define mem(a, b) memset(a, b, sizeof (a))
using namespace std;
///浮点数----------------------------------------------------------------------------------------------------------------------------
const db eps = 1e-9;
const db pi = acos(-1);
int sign(db x){ if(fabs(x) < eps) return 0;  return x > 0 ? 1 : -1; }
int cmp(db k1, db k2){ return sign(k1-k2); }
int inmid(db k1, db k2, db k3){ return sign(k2-k1)*sign(k3-k1) <= 0; }//k1在k2、k3之间
int intersect(db l1,db r1,db l2,db r2){// 区间相交判定,区间1在区间2前
    if(l1>r1) swap(l1,r1); if(l2>r2) swap(l2,r2); return cmp(r1,l2)!=-1 && cmp(r2,l1)!=-1;
}
///点、向量----------------------------------------------------------------------------------------------------------------------------
struct point{
    db x, y;
    point(){}
    point(db k1, db k2){ x = k1, y = k2; }
    point operator + (const point &k1) const { return point(x+k1.x, y+k1.y); }//向量加法、点+向量=点
    point operator - (const point &k1) const { return point(x-k1.x, y-k1.y); }//向量减法、点-点=向量
    point operator * (db k1) const { return (point){x*k1, y*k1}; }//向量数乘
    point operator / (db k1) const { return (point){x/k1, y/k1}; }//向量数除
    int operator == (const point &k1) const { return cmp(x,k1.x)==0 && cmp(y,k1.y)==0; }//比较两个点(向量)是否相同
    point turn(db k1){return (point){x*cos(k1)-y*sin(k1), x*sin(k1)+y*cos(k1)};}//逆时针旋转
    point turn90(){return (point){-y, x};}//逆时针旋转90度
    bool operator < (const point k1) const{
        int a = cmp(x, k1.x);
        if(a == -1) return 1; else if(a == 1) return 0; else return cmp(y,k1.y)==-1;
    }//比较两个点(向量)的大小,x越小则点越小;若x相等,则y越小点越小。可以实现按点的坐标排序
    db len(){ return sqrt(x*x+y*y); }//向量模长
    db len2(){ return x*x+y*y; }//向量模长的平方
    point dir(){ return (*this)/(*this).len(); }
    point unit(){ db k = len(); return point(x/k, y/k); }//单位向量
    db angle() { return atan2(y, x); }//向量的极角
    point getdel(){if (sign(x)==-1||(sign(x)==0&&sign(y)==-1)) return (*this)*(-1); else return (*this);}
    //当横坐标为负时,或横坐标为0纵坐标为负时,将点按原点做对称??? 角度是[-π/2, π/2]
    int getp() const { return sign(y)==1 || (sign(y)==0&&sign(x)==-1); }
    //判断点在12象限,或者在x的负半轴上??? 角度是(0, π]
    void scan(){ scanf("%lf%lf", &x, &y); };
    void print(){ printf("%.11f %.11f\n", x, y); }
};
int inmid(point k1, point k2, point k3){ return inmid(k1.x,k2.x,k3.x) && inmid(k1.y,k2.y,k3.y); }//k1 在 [k2,k3] 内
point midpo(point k1, point k2){ return (k1+k2)/2; }//得到两点中点
db dis2(point k1, point k2){ return (k1.x-k2.x)*(k1.x-k2.x) + (k1.y-k2.y)*(k1.y-k2.y); }
db dis(point k1, point k2){ return sqrt(dis2(k1, k2)); }
db cross(point k1, point k2){ return k1.x*k2.y - k1.y*k2.x; }//叉乘
db dot(point k1, point k2){ return k1.x*k2.x + k1.y*k2.y; }//点乘

struct line{// p[0]->p[1]
    point p[2];
    line(){}
    line(point k1,point k2){p[0]=k1; p[1]=k2;}
    point& operator [] (int k){return p[k];}
    int include(point k){ return sign(cross(p[1]-p[0],k-p[0]))>=0; }//点在直线左侧的判定
    point dir(){return p[1]-p[0];}//方向向量
    line push(){ // 向外 ( 左 ) 平移 eps
        point delta=(p[1]-p[0]).turn90().unit()*eps;
        return {p[0]-delta, p[1]-delta};
    }
};
point proj(point q, point k1, point k2){ point k=k2-k1; return k1+k*(dot(q-k1,k)/k.len2()); }// q 到直线 k1,k2 的投影
int checkonl(point q,point k1,point k2){ return sign(cross(k1-q, k2-k1))==0; }

struct circle{
    point o; db r;
    circle(){}
    circle(point _o, db _r){ o = _o, r = _r; }
    bool include(point k){ return cmp(dis(o, k), r) <= 0; }//点在圆内判定
};
vector<point> getcl(circle k1,point k2,point k3){ //求直线与圆的交点 沿着 k2->k3 方向给出 , 相切给出两个
    point k=proj(k2,k3,k1.o);
    db d=k1.r*k1.r-(k-k1.o).len2();
    if (sign(d)==-1) return {};
    point del=(k3-k2).unit()*sqrt(max((db)0.0, d));
    return {k-del, k+del};
}
int checkposcc(circle k1,circle k2){// 返回两个圆的公切线数量
    if (cmp(k1.r,k2.r)==-1) swap(k1,k2);
    db d=dis(k1.o,k2.o);  int w1=cmp(d,k1.r+k2.r),w2=cmp(d,k1.r-k2.r);
    if (w1>0) return 4;//相离
    else if (w1==0) return 3;//相切
    else if (w2>0) return 2; //相交
    else if (w2==0) return 1; //内切
    else return 0;//内含
}
vector<point> getcc(circle k1,circle k2){//求两圆交点 沿圆 k1 逆时针给出 , 相切给出两个
    int pd=checkposcc(k1,k2);
    if(pd==0||pd==4) return {};
    db a=(k2.o-k1.o).len2();
    db cosA=(k1.r*k1.r+a-k2.r*k2.r)/(2*k1.r*sqrt(max(a,(db)0.0)));
    db b=k1.r*cosA;
    db c=sqrt(max((db)0.0,k1.r*k1.r-b*b));
    point k=(k2.o-k1.o).unit(), m=k1.o+k*b, del=k.turn90()*c;
    return {m-del, m+del};
}
vector<point> tangentcp(circle k1,point k2){// 过圆外一点作圆的切线,求切点 沿圆 k1 逆时针给出
    db a=(k2-k1.o).len(),b=k1.r*k1.r/a,c=sqrt(max((db)0.0,k1.r*k1.r-b*b));
    point k=(k2-k1.o).unit(),m=k1.o+k*b,del=k.turn90()*c;
    return {m-del, m+del};
}
vector<line> tangentoutcc(circle k1,circle k2){//求两圆的外切线
    int pd=checkposcc(k1,k2); if (pd==0) return {};
    if (pd==1){//内含,返回一条切线
        point p1=getcc(k1,k2)[0]; point p2=p1+((p1-k1.o).turn90()/(p1-k1.o).len());
        return {(line){p1,p2}};
    }
    if (cmp(k1.r,k2.r)==0){
        point del=(k2.o-k1.o).unit().turn90().getdel();
        return {(line){k1.o-del*k1.r,k2.o-del*k2.r},(line){k1.o+del*k1.r,k2.o+del*k2.r}};
    } else {
        point p=(k2.o*k1.r-k1.o*k2.r)/(k1.r-k2.r);
        vector<point>A=tangentcp(k1,p),B=tangentcp(k2,p);
        vector<line>ans; for (int i=0;i<A.size();i++) ans.push_back((line){A[i],B[i]});
        return ans;
    }
}
vector<line> tangentincc(circle k1,circle k2){//求两圆的内切线
    int pd=checkposcc(k1,k2); if (pd<=2) return {};
    if (pd==3){
        point p1=getcc(k1,k2)[0]; point p2=p1+((p1-k1.o).turn90()/(p1-k1.o).len());
        return {(line){p1,p2}};
    }
    point p=(k2.o*k1.r+k1.o*k2.r)/(k1.r+k2.r);
    vector<point>A=tangentcp(k1,p),B=tangentcp(k2,p);
    vector<line>ans; for (int i=0;i<A.size();i++) ans.push_back((line){A[i],B[i]});
    return ans;
}
vector<line> tangentcc(circle k1,circle k2){//求两圆所有切线
    int flag=0; if (k1.r<k2.r) swap(k1,k2),flag=1;
    vector<line>A=tangentoutcc(k1,k2),B=tangentincc(k1,k2);
    for (line k:B) A.push_back(k);
    if (flag) for (line &k:A) swap(k[0],k[1]);
    return A;
}

circle getinvertcir(circle c1, circle c0){//c1关于c0的反演圆
    circle c2;
    db x0 = c0.o.x, y0 = c0.o.y, r0 = c0.r, x1 = c1.o.x, y1 = c1.o.y, r1 = c1.r;
    db d01 = dis(c0.o, c1.o);
    c2.r = 0.5*((1/(d01-r1))-(1/(d01+r1)))*r0*r0;
    db d02 = r0*r0/(d01+r1)+c2.r;
    //db _d02 = r0*r0/(d01-r1)-c2.r;
    c2.o.x = x0 + d02/d01*(x1-x0);
    c2.o.y = y0 + d02/d01*(y1-y0);
    return c2;
}

circle getinvertcir(line k1, circle c0){//直线k1关于c0的反演圆
    point a = proj(c0.o, k1[0], k1[1]);
    db oa = dis(c0.o, a);
    db ob = c0.r*c0.r/oa;
    point v = a-c0.o;  v = v/v.len();
    point b = c0.o+v*ob;
    circle res;
    res.o = midpo(c0.o, b);
    res.r = ob/2;
    return res;
}

line getinvertline(circle c1, circle c0){//c1关于c0的反演直线
    point v = c1.o - c0.o;
    v = v / v.len();
    db d = c0.r*c0.r / (2*c1.r);
    point k1 = v * d + c0.o;
    v = v.turn90();
    point k2 = k1 + v;
    return line(k1, k2);
}

point getinvertpoint(point p, circle c){//p关于c的反演点
    point v = (p-c.o).dir();
    db len = c.r*c.r/dis(c.o, p);
    return c.o+v*len;
}

circle c[10];

bool check(circle k){
    if(checkposcc(k, c[1]) != 3) return false;
    if(checkposcc(k, c[2]) != 3) return false;
    if(cmp(dis(k.o, c[0].o), k.r) != 0) return false;
    return true;
}

int main(){
    int T;
    cin >> T;
    while(T--){
        cin >> c[1].o.x >> c[1].o.y >> c[1].r;
        cin >> c[2].o.x >> c[2].o.y >> c[2].r;
        cin >> c[0].o.x >> c[0].o.y;
        c[0].r = 15;
        c[3] = getinvertcir(c[1], c[0]);
        c[4] = getinvertcir(c[2], c[0]);
        vector<line>A = tangentcc(c[3], c[4]);
        vector<circle>ans;
        for(int i = 0; i < A.size(); i++){
            if(checkonl(c[0].o, A[i][0], A[i][1])) continue;
            circle tmp = getinvertcir(A[i], c[0]);
            if(check(tmp)) ans.pb(tmp);
        }
        cout << ans.size() << endl;
        for(int i = 0; i < ans.size(); i++)
            printf("%.8f %.8f %.8f\n", ans[i].o.x, ans[i].o.y, ans[i].r);
    }
    return 0;
}

 

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值