【HDU 4773】Problem of Apollonius(圆的反演)

传送门

定一个常数 R R R和反演中心 O O O
反演就是对于每一个点 A A A变换到 A ′ A' A满足 O A ∗ O A ′ = R 2 OA*OA'=R^2 OAOA=R2,其中 A , A ′ , O A,A',O A,A,O共线
对于反演有如下性质:
1、圆反演之后还是一个圆(过反演中心的变成一条直线,视为特殊的圆)
2、一条直线反演之后为过 O O O的圆
3、反演之后圆之间关系(相交 / 切 / 离)不变

圆的反演
设反演中心为 O O O,圆心为 C 1 C_1 C1,半径为 r 1 r_1 r1,常数为 r r r
那么反演之后的圆心为 C 2 C_2 C2,半径为 r 2 r_2 r2

图是网上搬的,侵删

那么有
( O C 1 + r 1 ) ( O C 2 − r 2 ) = r 2 ( O C 1 − r 1 ) ( O C 2 + r 2 ) = r 2 (OC_1+r_1)(OC_2-r_2)=r^2\\ (OC_1-r1)(OC_2+r2)=r^2 (OC1+r1)(OC2r2)=r2(OC1r1)(OC2+r2)=r2
那么 r 2 = 1 2 r 2 ( 1 O C 1 − r 1 − 1 O C 1 + r 1 ) r_2=\frac 1 2r^2(\frac{1}{OC_1-r_1}-\frac 1{OC_1+r_1}) r2=21r2(OC1r11OC1+r11)
O C 2 = O A ′ + r 2 = 1 2 r 2 ( 1 O C 1 − r 1 + 1 O C 1 + r 1 ) OC_2=OA'+r_2=\frac 1 2r^2(\frac{1}{OC_1-r_1}+\frac 1{OC_1+r_1}) OC2=OA+r2=21r2(OC1r11+OC1+r11)
C 2 = O + ( C 1 − O ) ∗ ( O C 2 / O C 1 ) C_2=O+(C_1-O)*(OC_2/OC_1) C2=O+(C1O)(OC2/OC1)


题意:给定两个圆 A , B A,B A,B和一个点 O O O,求所有满足过 O O O且与两圆外切的圆

O O O为反演中心反演
那么就变成反演之后两个圆的外公切线
O O O与两个圆在同一侧
求公切线可以先求出公切线的倾角

然后考虑直线转回圆
显然 2 r ∗ d i s ( 反 演 中 心 到 直 线 ) = R 2 2r*dis(反演中心到直线)=R^2 2rdis(线)=R2
圆心可以用反演中心到圆心的向量与 A A A和切点的向量相似求出

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define pb push_back
#define pii pair<int,int>
#define ll long long
#define bg begin
#define fi first
#define se second
cs int RLEN=1<<20|1;
inline char gc(){
    static char ibuf[RLEN],*ib,*ob;
    (ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
    return (ib==ob)?EOF:*ib++;
}
inline int read(){
    char ch=gc();
    int res=0;bool f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
    return f?res:-res;
}
inline ll readll(){
    char ch=gc();
    ll res=0;bool f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
    return f?res:-res;
}
inline int readstring(char *s){
    int top=0;char ch=gc();
    while(isspace(ch))ch=gc();
    while(!isspace(ch)&&ch!=EOF)s[++top]=ch,ch=gc();
    return top;
}
template<typename tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<typename tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
namespace Module{
cs int mod=998244353;
inline int add(int a,int b){return (a+b>=mod)?(a+b-mod):(a+b);}
inline int dec(int a,int b){return (a<b)?(a-b+mod):(a-b);}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return (r>=mod)?(r%mod):r;}
inline void Add(int &a,int b){a=(a+b>=mod)?(a+b-mod):(a+b);}
inline void Dec(int &a,int b){a=(a<b)?(a-b+mod):(a-b);}
inline void Mul(int a,int b){static ll r;r=1ll*a*b;a=(r>=mod)?(r%mod):r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
inline int fix(ll x){x%=mod,(x<0)&&(x+=mod);return x;}
}
cs int N=200010;
struct pt{
    double x,y;
    pt(double _x=0,double _y=0):x(_x),y(_y){}
    friend inline pt operator +(cs pt &a,cs pt &b){
        return pt(a.x+b.x,a.y+b.y);
    }
    friend inline pt operator -(cs pt &a,cs pt &b){
        return pt(a.x-b.x,a.y-b.y);
    }
    friend inline double operator *(cs pt &a,cs pt &b){
        return a.x*b.y-a.y*b.x;
    }
    friend inline double operator ^(cs pt &a,cs pt &b){
        return a.x*b.x-a.y*b.y;
    }
    friend inline pt operator *(cs pt &a,cs double &b){
        return pt(a.x*b,a.y*b);
    }
    friend inline pt operator /(cs pt &a,cs double &b){
        return pt(a.x/b,a.y/b);
    }
    inline pt move(double a,double d){return pt(x+d*cos(a),y+d*sin(a));}
    inline double ang(){return atan2(y,x);}
    inline double dis(){return sqrt(x*x+y*y);}
}P;
cs double eps=1e-10,INF=1e12;
inline int sign(double x){
    return (x>eps)-(x<-eps);
}
struct cir{
    pt o;double r;
    cir(){}
    cir(pt x,double _r){o=x,r=_r;}
    inline void write(){
        printf("%.8lf %.8lf %.8lf\n",o.x,o.y,r);
    }
}A,B,ans[4];
cs double R=6;
int tot;
inline cir CInv(cir A){
    cir res;
    double oc1=(A.o-P).dis();
    double v1=1.0/(oc1-A.r),v2=1.0/(oc1+A.r),oc2=0.5*R*R*(v1+v2);
    res.r=0.5*R*R*(v1-v2),res.o=P+(A.o-P)*(oc2/oc1);
    return res;
}
inline void getans(pt p1,pt p2){
    tot++;
    double di=fabs((P-p1)*(P-p2))/(p1-p2).dis()*2.0;
    ans[tot].r=R*R/di;
    double d=(p1-A.o).dis();
    ans[tot].o=P+(p1-A.o)*(ans[tot].r/d);
}
inline void solve(){
    tot=0;
    A=CInv(A),B=CInv(B);
    pt del=B.o-A.o;
    double a1=del.ang(),a2=acos((A.r-B.r)/(A.o-B.o).dis());
    pt p1=A.o.move(a1+a2,A.r),p2=B.o.move(a1+a2,B.r);
    if(sign((p2-p1)*(A.o-p1))==sign((p2-p1)*(P-p1)))getans(p1,p2);
    p1=A.o.move(a1-a2,A.r),p2=B.o.move(a1-a2,B.r);
    if(sign((p2-p1)*(A.o-p1))==sign((p2-p1)*(P-p1)))getans(p1,p2);
}
int main(){
    #ifdef Stargazer
    freopen("lx.in","r",stdin);
    #endif
    int T=read();
    while(T--){
        A.o.x=read(),A.o.y=read(),A.r=read();
        B.o.x=read(),B.o.y=read(),B.r=read();
        P.x=read(),P.y=read();
        solve();
        cout<<tot<<'\n';
        for(int i=1;i<=tot;i++)ans[i].write();
    }return 0;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值