前言

今天的python课，我给学生介绍math库。其中一个孩子问，能看到PI的值吗？我马上开始了我的表演。进入命令行，导入，调用math.pi,一气呵成。谁知，这孩子并不买账！嚷嚷道：我要看小数点后100位@@@

一、马青公式

Machin公式：π=16arctan(1/5)-4arctan(1/239)

二、Python代码实现

1.理清思路

计算π 的公式已经拿到，我们有一个明确的想要精确的π的位数，那么当式子中n越大的时候，第n项就会越小，小到计算这个精确目标时可以忽略的地步。

注意1：由于程序中浮点数的有效位数有限（小数点后很长的话会被省略，很大的数如果是小数也会用科学技术法），会损失很大精度，所以用整数来计算π。

比如想计算小数点后 7 位，就计算 π * 10^7 是多少，返回 31415926。比如想计算小数点后 3 位，程序就会给你返回整数 3141，是扩大了 10^3 的结果。

注意2: 想用整数计算 π，那除法的时候就要使用取模运算来保证结果也是整数，以达到保证精度的目的。

注意3: 用取模操作代替除法，但是取模对于除法来说本身就是有精度损失的，为了弥补这部分精度，可以先把想精确的位数增加 k，最后把结果减少 k 位。

注意4: 注意n:求 π 的数学公式中的 n，还是精确到π 的位数。

代码如下（python）：

# 计算准确圆周率的马青公式
import time
n = int(input('请输入圆周率小数点后的位数:'))
start_time = time.time()   # 计算pi计时起点
w = n + 10                 # 考虑精度，位数增加 k，结果减少 k 位
b = 10 ** w                # 用整数计算pi  如3.14 * 10^2 = 314
x1 = b * 4// 5             # 第一项前半部分
x2 = b // -239             # 第一项后半部分
sum = x1 + x2              # 第一项的值
n *= 2
for i in range(3, n, 2):x1 //= -25x2 //= -239*239x = (x1+x2)//isum += x
mpi = sum * 4
mpi //= 10**10
end_time = time.time()
run_time = str(end_time - start_time)
mpi_str = str(mpi)
mpi_str = mpi_str[:1] + '.' + mpi_str[1:]
# f = open('mpi.txt', 'w')
# f.write(mpi_str)
# f.close()
print('运行时间: ' + run_time)
print('计算结果: ',mpi_str)
print('')

2.C++版本

代码如下（C++）：

/***********************************************************使用梅钦公式π=16arctan(1/5)-4arctan(1/239)*通过建立数组来储存超过类型范围的小数部分*用long类型数组pi[ ]来储存各项（每一项有4位)*使用long类型term[ ]来实现对每一次泰勒展开结果的储存*之后通过进位等操作，实现用数组储存足够位数的π。         ***********************************************************/
void machin_pi(long L)
{//常量C用于进位退位操作，10000代表pi数组以四位进行储存const long C = 10000;//flag用于正负号的判定short flag;//梅钦公式具体应用时π=80*（1/25-1/（25*25*3）+......)+(-956)*(......)//xishu[]为两项展开的外系数，div[]为每次进阶指数所用的除数//odd为1,3，5等奇数，shang_1，r_1是进阶指数的商和余数,r_2是计算π各项的余数long xishu[2] = { 80,-956 }, div[2] = { 25,57121 },odd,shang_1, r_1, r_2,i=0,j,t,tt,k=0,len=L;L = L / 4 + 3;//term[]用来完成对每个升幂项的足够长度保存//pi[]用来储存最后结果long *term = new long[L];long *pi = new long[L];while (i < L)pi[i++] = 0;while (k - 2){//完成对term的初始化flag = 1;term[0] = xishu[k];i = 1;while (i < L){term[i] = 0;++i;}for (i = 0,odd=1; i < L; ++i){j = i;for (r_1=r_2 = 0;j < L; ++j){t = r_1 * C + term[j];                //t、r_1、term[]用来储存25与57121的各幂term[j] = t / div[k];                 //每次进阶都用上一轮已经算好的term数组进行升阶（后一项都是前一项除以25或57121）r_1 = t % div[k];                     //保存余数，在下一项中×C用于退位tt = r_2 * C + term[j];               //tt、r_2、shang_1用来储存泰勒展开各项shang_1 = tt / odd;                   //在term[]的基础上除以奇数1、3、5等r_2 = tt % odd;                       //保存余数，在下一项中用于退位pi[j] += (flag?shang_1:-shang_1);     //用flag进行标定正负号}if (term[i])--i;                          //当term[i]==0时，后续展开不会影响这一项（这四位小数），进行++i,直接从下一项开始odd += 2; flag=!flag;                               //泰勒展开正负号改变                               }++k;                                          //处理完arccot5后进入下一轮处理arccot239}k = 0;//上述循环之后，pi[]中已经储存了各项数据//但pi[]中数据可能存在负数、五位数及以上等等情况//接下来的循环进行处理，通过进位全部转化成正四位数while (--i){pi[i] += k;if (pi[i] < 0){k = pi[i] / C-1;pi[i] %= C;pi[i] += C;}else{k = pi[i] / C;pi[i] %= C;}}//输出圆周率printf("3.");for (i = 1; i < len / 4 + 1; ++i){printf("%04ld", pi[i]);}switch (len%4){case 0:break;case 1:printf("%01ld", pi[i]/ 1000); break;case 2:printf("%02ld", pi[i]/ 100); break;case 3:printf("%03ld", pi[i]/ 10); break;}//释放内存delete[] term;delete[] pi;
}

总结

注意3中提出的这个方案是有漏洞的，因为 k 是被固定在程序中的，随着用户输入 n 的增大，越来越大量的取模运算损失的精度，必定是固定不变的 k 值所不能弥补的。所以对每一个固定的 k，一定会有一个 n 值，使得本程序所计算的结果是错误的。

对此，我还有一个方案，k 值不固定，由 n 计算得出，n 越大就取越大的 k。但我不想随便定一个 k 与 n 的关系，当然通过代码可以试出一个不错的，但程序计算总有尽头，没有理论支持就无法证明这个关系是真正合理的。当 k 取 10 时，n 取 100000 还都没有出现误差。

【Python】计算任意位数的圆周率π（Machin Formula）相关推荐

【Python、数学】计算任意位数的圆周率π（马青公式）
1. 公式准备计算准确圆周率的马青公式: 对反正切进行级数展开: 就可以得到 π = 16(1/5 - 1/3/5^3 + 1/5/5^5 - ...) - 4(1/239 - 1/3/239^3 ...
python计算圆周率近似值_使用MicroPython计算任意位数圆周率
计算任意精度的圆周率是个有趣的主题,得益于python的强大计算能力,我们在MicroPython中也可以轻松的计算pi的数值. 先输入下面的代码: """ 文件:pi. ...
python计算任意多边形面积
看了一些方法网上的资料,也自己倒腾了一些想法. 然后看到这个大神的博客(<计算任意多边形的面积>),我就服气了. 我把代码转换为python: # 计算任意多边形的面积,顶点按照顺时针或者 ...
Python 计算任意两向量之间的夹角
如图所示,我们要计算任意两个向量之间的夹角. (图中的坐标数字是估计值,随手给定) python代码如下 import math AB = [1,-3,5,-1] CD = [4,1,4.5,4.5] ...
python 计算数字位数,Python | 计算一个数字的总位数
先决条件: Python中的二进制数系统给定一个数字,我们必须使用Python查找二进制值的总位数来表示该数字. 示例 Input: num = 61 Binary value of 61 is = ...
Python计算任意数据的分布函数(累积概率函数CDF)与概率密度(PDF)，并绘图
分布函数求导就是概率密度,在任意的数据中由于数据分布不确定,所以严格来讲这些数据是没有分布函数和概率密度函数的,但是我们就是想得出那种结果,画出那种曲线,这里就需要借助scipy.stats.relf ...
python 计算数字位数_Python 统计位数为偶数的数字的代码
问题描述给出一个整数数组 nums,请返回其中位数为偶数的数字的个数. 示例 1: 输入:nums = [12,345,2,6,7896] 输出:2 解释: 12 是 2 位数字(位数为偶数) 34 ...
Python计算任意单调曲线在给定区间上的近似长度
本文要点在于:任意曲线在一个很小的局部都可以看作直线. def curveLength(xs, func): '''xs:x轴的采样点,越密越准确 func:曲线方程对应的函数''' #函数曲线上的采 ...
Python计算开方、立方、圆周率，精确到小数点后任意位
Python计算的位数在电脑上做了一个实验,看看python能计算到多少位,一下是结果. x = math.sqrt((3)) print ("%.53f"%(x)) print ...

【Python】计算任意位数的圆周率π（Machin Formula）

前言

一、马青公式

二、Python代码实现

1.理清思路

2.C++版本

总结

【Python】计算任意位数的圆周率π（Machin Formula）相关推荐

最新文章

热门文章