查看原文
其他

利用多边形估计Pi值

Y叔叔 YuLabSMU 2022-09-20


这是想翘动地球的阿基米德所提出的方法,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 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

往期精彩

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存