题面
这题有两种方法(然而两种我都想不到)
方法一
前置芝士
笛卡尔定理
我们定义一个圆的曲率为(k=pm {1over r}),其中(r)是圆的半径
若在平面上有两两相切,且六个切点互不相同的四个圆,设其曲率分别为(k1,k2,k3,k4)(若该圆和其它所有圆都外切,则其曲率取正,否则曲率取负),则有
类似的,若是空间中有两两相切且切点互不相同的五个球体,则有
题解
首先那个粉色圆的直径就是(R-r),于是我们现在知道前三个圆的半径了
上面的柿子可以看做一个关于(k4)的方程,化简一下可得
这方程一看就很不好解的样子……
想想我们初中时候学过的韦达定理,对于方程(ax^2+bx+c=0),设它的两个解为(x_1,x_2),则有
设(k3)代表绿圆(i)的曲率,(k2)代表黄圆的曲率,(k1)代表最大的圆的曲率,我们要求的就是与这三个圆都相切的圆的曲率……
但是这里是两个解的和的形式怎么办丫……
冷静想想,和这三个圆都相切的圆,除了绿圆(i+1),似乎只剩下绿圆(i-1)?而(i-1)的曲率我们之前已经知道了?
那么我们就可以(O(n))递推了,时间复杂度(O(nT))
//minamoto
#include<bits/stdc++.h>
#define R register
#define inline __inline__ __attribute__((always_inline))
#define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
R int res;R char ch;
while((ch=getc())>'9'||ch<'0');
for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
return res;
}
int n,R1,R2,R3;double r1,r2,r3,r4,r5;
int main(){
// freopen("testdata.in","r",stdin);
for(int T=read();T;--T){
R1=read(),R2=read(),n=read(),R3=R1-R2;
r1=-1.0/R1,r2=1.0/R2,r3=1.0/R3,r4=r1+r2+r3;
fp(i,2,n)r5=r3,r3=r4,r4=(r1+r2+r3)*2-r5;
printf("%.10lf
",1.0/r4);
}
return 0;
}
方法二
前置芝士
圆的反演
题解
我们设黄圆为(A),绿圆为(B),大圆就叫大圆,并且以下假设(A)和大圆的交点为原点,这三个圆共同的直径为(x)轴
如果我们以原点为反演中心,那么反演之后(A)和大圆会变成两条直线。又因为反演不改变相切关系,所以(B)反演之后就在这两条直线中间。显然(B)反演之后的圆心仍在(x)轴上
如果我们还要放上别的圆(假设一直往上放),由于它与这三个圆有唯一交点,所以肯定是被卡在两条直线中间,且刚好在之前的圆上方
因为这些圆的半径相同,我们可以(O(1))得到第(n)个圆的圆心的坐标。之后我们连接圆心和原点,这条直线和圆的两个交点分别记为(p1,p2),把它们反演回来,对应的两个点之间的距离就是原来的第(n)个圆的半径了
//minamoto
#include<bits/stdc++.h>
#define R register
#define inline __inline__ __attribute__((always_inline))
#define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
R int res,f=1;R char ch;
while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
return res*f;
}
struct node{
double x,y;
inline node(){}
inline node(R double xx,R double yy):x(xx),y(yy){}
inline node operator +(const node &b)const{return node(x+b.x,y+b.y);}
inline node operator -(const node &b)const{return node(x-b.x,y-b.y);}
inline node operator *(const double &b)const{return node(x*b,y*b);}
inline double norm(){return sqrt(x*x+y*y);}
}p,pp,o;
int n,r1,r2;double r,lx,rx,ir;
void get(R node o,R double k,R double b){
double del=ir/sqrt(k*k+1);
p=node(o.x-del,k*(o.x-del)+b),
pp=node(o.x+del,k*(o.x+del)+b);
}
inline node inv(R node p){
double len=1.0/p.norm();
return p*(r*r*len*len);
}
int main(){
// freopen("testdata.in","r",stdin);
for(int T=read();T;--T){
r1=read(),r2=read(),n=read(),r=r2*0.5;
lx=r*r/(r1*2),rx=r*r/(r2*2),ir=(rx-lx)*0.5;
o=node((lx+rx)*0.5,2*n*ir);
get(o,o.y/o.x,0);
printf("%.10lf
",(inv(p)-inv(pp)).norm()*0.5);
}
return 0;
}