题面
这题有两种方法（然而两种我都想不到）
方法一
前置芝士
笛卡尔定理
我们定义一个圆的曲率为(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;
}