传送门:http://acm.hdu.edu.cn/showproblem.php?pid=6158
本题是一个计算几何题——四圆相切。
平面上的一对内切圆,半径分别为R和r。现在这一对内切圆之间,按照如图所示的方式依次放置N个相切的圆。求放置的这N个圆的面积之和。
在此,首先介绍一个定理:笛卡尔定理。Wiki: Descartes' theorem。
平面上的四个圆,第i个圆的半径为r[i],曲率为κ[i](注:κ=r-1)。若这四个圆中的每一对均构成外切,则其曲率满足约束:
$\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-1)。半径的约束式相应地写成:
$\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(-1); const double eps = 1e-14; 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", 0); 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(0); double ans = rad * rad; //Add following circles. for (int i = 2; i <= n; i += 2) { double rad = 1.0 / get_curv(i / 2); double ds = rad * rad; if (ds < eps) break; ans += ds * (i < n ? 2 : 1); } printf("%.5f\n", ans * pi); } return 0; }
本题还有一种更为简单的解法,即通过笛卡尔定理与韦达定理进行迭代。参考程序如下:
#include <bits/stdc++.h> using namespace std; const double pi = acos(-1); const double eps = 1e-14; 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", 0); 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 = 2; i <= n; i += 2) { double ds = 1.0 / (k_4 * k_4); if (ds < eps) break; ans += ds * (i < n ? 2 : 1); 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 0; }