[C++]利用多边形估计Pi值
这是想翘动地球的阿基米德所提出的方法，O圆的中心，R是半径，先从6边形开始，再切12边形，24边形，48边形，不断切之后，多边形的周长会越来越接近周长，就可以拿来估计Pi值。
比如上面这个图，蓝色这条线长度为D，是多边线的边，现在在它中间一切为二，现在我们要求新的边x的长度。蓝色的一半是d=D/2, 新的顶点到蓝色线是y,里面那段是h = R-y。
那么很容易推出：
h = sqrt(R^2 - d^2)
x = sqrt(d^2 + (R-h)^2)
所以呢，这里x是d的函数，而6边形D=R，我们从6边形开始，不断对切，算出新的x出来，就可以拿来算Pi值。
因为小数点计算常有精度问题，我们这里使用大整数来算，假设R是1e18，内接个96边形，侧边长是65438165643552283，那么Pi值就可以估计为：
Pi = 96 * 65438165643552283 / 2 = 3141031950890509584这种题最好用python或java来解，因为内置支持大整数。
这道题的输入有两个数，K和N，K用来指定R = 10^K，而N用来指定多边形的边数是6*2^N。
要用C++来解的话，需要额外的库，这里我会用GMP库。
#include <iostream>
#include <gmpxx.h>
struct polygon_t {
mpz_class R;
mpz_class d;
mpz_class h;
mpz_class side;
long int n;
};
mpz_class mpz_pow(mpz_class x, int n);
polygon_t dodecagon(mpz_class R);
polygon_t update_polygon(polygon_t poly);
polygon_t polygon(mpz_class R, int N);
int main() {
int K=57;
int N=30;
mpz_class R(1);
for (int i=0; i<K; i++) {
R *= 10;
}
polygon_t poly = polygon(R, N);
mpz_class Pi = poly.n * poly.side/2;
std::cout << Pi << std::endl;
}
先定义个结构体，放多边形的信息，主函数照常会简单，polygon函数传入R和N，返回多边形的信息，计算Pi，打印出来。
6边形的时候对切d=R/2，我们可以很容易算出来12边形，以此为基础。
polygon_t dodecagon(mpz_class R) {
polygon_t res;
res.R = R;
res.d = R/2;
res.h = sqrt(mpz_pow(R,2) - mpz_pow(res.d, 2));
res.side = sqrt(mpz_pow(res.d, 2) + mpz_pow(res.R-res.h, 2));
res.n = 12;
return res;
}
有了一个多边形之后，不断对切，需要通过上一个多边形，计算对切后的下一个多边形：
polygon_t update_polygon(polygon_t poly) {
polygon_t res;
res.R = poly.R;
res.d = poly.side/2;
res.h = sqrt(mpz_pow(res.R,2) - mpz_pow(res.d, 2));
res.side = sqrt(mpz_pow(res.d, 2) + mpz_pow(res.R-res.h, 2));
res.n = poly.n * 2;
return res;
}
不管是初始化的12边形，还是更新多边形，都是利用上面的x和d的关系，函数很简单。但是有了这两个函数之后，初始化一个12边形，不断对切，我们可以得到指定N的任意多边形的信息：
polygon_t polygon(mpz_class R, int N) {
polygon_t poly = dodecagon(R);
for (int i=1; i<N; i++) {
poly = update_polygon(poly);
}
return poly;
}
有了这个函数，就有了主程序的简单思路，来个多边形，估个Pi值。
PS 里面用到的mpz_pow是自己定义的：
mpz_class mpz_pow(mpz_class x, int n) {
mpz_class res(1);
for (int i=0; i<n; i++) {
res *= x;
}
return(res);
}
这道题全部放在主函数里，也不长。但通过divide-and-conquer，简单清晰了不少。
用上面K=57, N=30，算出来是：
3141592653589793238338135707214807701028989418292419493888
PS2 上面的程序需要用下面的指令来编译：
g++ -I/usr/local/opt/gmp/include -L/usr/local/opt/gmp/lib -lgmpxx -lgmp pi.cpp