传送门：http://acm.hdu.edu.cn/showproblem.php?pid=6158
本题是一个计算几何题——四圆相切。
平面上的一对内切圆，半径分别为R和r。现在这一对内切圆之间，按照如图所示的方式依次放置N个相切的圆。求放置的这N个圆的面积之和。
在此，首先介绍一个定理：笛卡尔定理。Wiki: Descartes' theorem。
平面上的四个圆，第i个圆的半径为r[i]，曲率为κ[i]（注：κ=r）。若这四个圆中的每一对均构成外切，则其曲率满足约束：
$\left(\sum_{i=1}^4 \kappa _i\right)^2 = 2\cdot \sum_{i=1}^4 \kappa _i^2$
通过这个定理，可以得到以下情景相应的约束：
平面上的三个圆，第i个圆的半径为r[i]，曲率为κ[i]。若这三个圆中的每一对均构成外切，且同时内切于一个半径为R，曲率为K的大圆，则其曲率同样满足以上的关系（注：此处大圆的曲率应取负值，即K=-R）。半径的约束式相应地写成：
$\left(\sum_{i=1}^3 \kappa _i -\frac{1}{R} \right)^2 = 2\left( \sum_{i=1}^3 \kappa _i^2+\frac{1}{R^2}\right)$
接下来，首先考虑上半侧的情况（下半侧与之对称）。设上半侧放置的第k个圆的曲率为c[k]（约定放置于中间的圆的曲率为c[0]），则其与半径为r的圆、放置的第k-1个圆相外切，并同时内切于半径为R的圆。根据四圆相切的关系写出约束式：
$\left(\frac{1}{r}-\frac{1}{R} +c_k+c_{k-1}\right)^2 = 2\left(\frac{1}{r^2}+\frac{1}{R^2} +c_k^2+c_{k-1}^2\right)$
相应地考虑第k+1个圆，则有：
$\left(\frac{1}{r}-\frac{1}{R} +c_k+c_{k+1}\right)^2 = 2\left(\frac{1}{r^2}+\frac{1}{R^2} +c_k^2+c_{k+1}^2\right)$
两式相减，则有：
$(c_{k+1}-c_{k-1})\left( 2\frac{R-r}{Rr}+2c_k+c_{k+1}+c_{k-1}\right )=2(c_{k+1}+c_{k-1})(c_{k+1}-c_{k-1})\\\Rightarrow 2\frac{R-r}{Rr}+2c_k=c_{k+1}+c_{k-1}\Rightarrow (c_{k+1}-c_k)-(c_k-c_{k-1})=2\frac{R-r}{Rr}$
设d[k]=c[k]-c[k-1]，则d[]是一个等差数列。为求得这个等差数列，首先需要求解首项。
c[0]是显然的，而c[1]则可以借助与R、r、c[0]的关系求解。
$c_0=\frac{1}{R-r}\\c_1=\frac{R^2+r^2-Rr}{Rr(R-r)}$
于是，d[]的通项公式：$d_k=(2k-1)\frac{R-r}{Rr},k=1,2,3,\cdots$
于是，c[]的通项公式：$c_k=\frac{1}{R-r}+\frac{R-r}{Rr}k^2,k=0,1,2,3,\cdots$
求解时注意精度控制。参考程序如下：
#include <bits/stdc++.h>
using namespace std; const double pi = acos(-);
const double eps = 1e-; int R, r;
double a, b; double get_curv(int k)
{
return a + b * k * k;
} int main(void)
{
int t;
scanf("%d", &t);
while (t--) {
int n;
scanf("%d%d%d", &R, &r, &n);
if (R == r) {
printf("%.5f\n", );
continue;
}
if (R < r) swap(R, r);
a = 1.0 / (R - r);
b = 1.0 * (R - r) / (R * r);
//Add first circle.
int rad = 1.0 / get_curv();
double ans = rad * rad;
//Add following circles.
for (int i = ; i <= n; i += ) {
double rad = 1.0 / get_curv(i / );
double ds = rad * rad;
if (ds < eps) break;
ans += ds * (i < n ? : );
}
printf("%.5f\n", ans * pi);
}
return ;
}
本题还有一种更为简单的解法，即通过笛卡尔定理与韦达定理进行迭代。参考程序如下：
#include <bits/stdc++.h>
using namespace std; const double pi = acos(-);
const double eps = 1e-; int R, r; int main(void)
{
int t;
scanf("%d", &t);
while (t--) {
int n;
scanf("%d%d%d", &R, &r, &n);
if (R == r) {
printf("%.5f\n", );
continue;
}
if (R < r) swap(R, r);
//Add first circle.
double ans = (R - r) * (R - r);
double k_1 = -1.0 / R;
double k_2 = 1.0 / r;
double k_3 = 1.0 / (R - r);
double k_4 = k_1 + k_2 + k_3;
//Add following circles.
for (int i = ; i <= n; i += ) {
double ds = 1.0 / (k_4 * k_4);
if (ds < eps) break;
ans += ds * (i < n ? : );
double k_5 = 2.0 * (k_1 + k_2 + k_4) - k_3;
k_3 = k_4;
k_4 = k_5;
}
printf("%.5f\n", ans * pi);
}
return ;
}