这是想翘动地球的阿基米德所提出的方法,O圆的中心,R是半径,先从6边形开始,再切12边形,24边形,48边形,不断切之后,多边形的周长会越来越接近周长,就可以拿来估计Pi值。

比如上面这个图,蓝色这条线长度为D,是多边线的边,现在在它中间一切为二,现在我们要求新的边x的长度。蓝色的一半是d=D/2, 新的顶点到蓝色线是y,里面那段是h = R-y。

那么很容易推出:

所以呢,这里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 xint n);  
polygon_t dodecagon(mpz_class R);  
polygon_t update_polygon(polygon_t poly);  
polygon_t polygon(mpz_class Rint 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