Machin formula /梅钦公式
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了Machin formula /梅钦公式相关的知识，希望对你有一定的参考价值。
梅钦公式( Machin formula )是一个计算 π 值的公式，至今仍被广泛应用，本文介绍如何使用计算机实现它。
π 的简史
四千年前，巴比伦人用 3+ 1/8 作为圆周率， 同时期的埃及人用 4-(8/9)^2 做圆周率；
三千年前，中国先人用3 作为圆周率；
公元前三世纪，古希腊科学家阿基米德首先采用计算的方法，得出 π 可能是 3.14；
公元五世纪，我国数学家祖冲之把 π 算到了 3.1415926 到 3.1415927 之间；
公元 15 世纪， 阿拉伯数学家阿尔·卡西 (Al-Kāshī,1380？– 1429)，用几何的方法，计算到了小数点后 16 位；
1666 年, 牛顿用自己设计的公式把 计算到了小数点后的15位,这个公式的收敛速度非常慢，我猜想他可能花了几个月，甚至几年干这事儿。
1706年，英国人 约翰·梅钦( John Machin) 发明了一个用于计算 π 值的公式。
1873 年, William Shanks 使用梅钦公式用了几年时间,计算到了 小数点后的 707 位。
20世纪30年代, 人们( 新西兰的数学家艾特肯(Aitken)教授 )才开始怀疑 Shanks 犯了一个错误。事实上,他在第528位犯了一个错误,所以剩下的180位数是错误的。
1949年，ENIAC（Electronic Numerical Integrator And Computer, 电子数字积分计算机）计算机用了70个小时,计算到了 小数点后的2037位。
近年来，被到了小数点后的1.24万亿位。
简介
虽然印度人拉马努金( Srinivasa Ramanujan ) 整了一堆如何计算 π 的公式，还有了 BBP ( Bailey- Borwein -Plouffe, David Bailey / Peter Borwein / Simon Plouffe )公式， 以及其他若干种类梅钦(Machin-like)公式，但梅钦公式至今仍然是计算 π 值的主要公式。
我已建立了百度词条 “梅钦公式”，请自行查看它是什么样的，因为：这里无法发数学公式。
梅钦公式是格里高利/莱布尼茨计算 公式的变体，但是更实用，它的收敛速度显著增加，这使得它成为了更实用的计算π的方法。
实现
多种编程语言可实现使用梅钦公式计算 π 小数点后任意精度的值，这里给出一种 javascript 实现：
/*
* @Author: coolwp.com
* @Date: 2017-08-22 05:40:48
* @Last Modified by: suifengtec
* @Last Modified time: 2017-08-23 03:24:11
**/
/*
用JavaScript 和 梅钦公式计算 PI。
使用
node pi.js
*/
var getPi = (function () {
function getPi(space) {
this.aX = [];
this.cellSize = 11;
this.tStart = new Date();
this.spaceStr = space ? space : " ";
this.aBigInt = Math.pow(10, this.cellSize);
}
getPi.prototype.getIt = function (length) {
var digitsLenAfterDot = Number(length);
var cellSize = this.cellSize;
var coeff = Array(10),
iAng = Array(10);
var arrLen = Math.ceil(1 + digitsLenAfterDot / cellSize);
var aPI = Array(arrLen);
var aArctan = Array(arrLen);
coeff[0] = 4;
coeff[1] = -1;
coeff[2] = 0;
iAng[0] = 5;
iAng[1] = 239;
iAng[2] = 0;
aPI = this.makeArr(arrLen, aPI, 0);
for (var i = 0; coeff[i] != 0; i++) {
aArctan = this.getArcTan(iAng[i], arrLen, aArctan);
aArctan = this.Mul(arrLen, aArctan, Math.abs(coeff[i]));
if (coeff[i] > 0) {
aPI = this.Add(arrLen, aPI, aArctan);
}
else {
aPI = this.Sub(arrLen, aPI, aArctan);
}
}
aPI = this.Mul(arrLen, aPI, 4);
var sPI = "3.";
var tempPI = "";
for (var i = 0; i < aPI.length; i++) {
aPI[i] = String(aPI[i]);
if (aPI[i].length < cellSize && i != 0) {
while (aPI[i].length < cellSize) {
aPI[i] = "0" + aPI[i];
}
}
tempPI += aPI[i];
}
for (var i = 0; i + 1 <= digitsLenAfterDot; i++) {
if (i == 0) {
sPI += tempPI.charAt(i + 1);
}
else {
var newLineSymbol = this.spaceStr, spForPer5Digits = this.spaceStr;
if ((i + 1) % 50 == 0 && i != 0) {
sPI += tempPI.charAt(i + 1) + newLineSymbol;
}
else {
if ((i + 1) % 5 == 0) {
sPI += tempPI.charAt(i + 1) + spForPer5Digits;
}
else {
sPI += tempPI.charAt(i + 1);
}
}
}
}
var tShutdown = new Date();
var timeTaken = (tShutdown.getTime() - this.tStart.getTime()) / 1000;
var timeSp = "耗时: " + timeTaken + " 秒";
return sPI + " \n" + timeSp;
};
getPi.prototype.Mul = function (n, aX, iMult) {
var carry = 0;
var prod;
for (var i = n - 1; i >= 0; i--) {
prod = (aX[i]) * iMult;
prod += carry;
if (prod >= this.aBigInt) {
carry = Math.floor(prod / this.aBigInt);
prod -= (carry * this.aBigInt);
}
else {
carry = 0;
}
aX[i] = prod;
}
return aX;
};
getPi.prototype.Div = function (n, aX, iDiv, aY) {
var carry = 0;
var currVal;
var theDiv;
for (var i = 0; i < n; i++) {
currVal = Number(aX[i]) + Number(carry * this.aBigInt);
theDiv = Math.floor(currVal / iDiv);
carry = currVal - theDiv * iDiv;
aY[i] = theDiv;
}
return aY;
};
getPi.prototype.Add = function (n, aX, aY) {
var carry = 0;
for (var i = n - 1; i >= 0; i--) {
aX[i] += Number(aY[i]) + Number(carry);
if (aX[i] < this.aBigInt) {
carry = 0;
}
else {
carry = 1;
aX[i] = Number(aX[i]) - Number(this.aBigInt);
}
}
return aX;
};
getPi.prototype.Sub = function (n, aX, aY) {
for (var i = n - 1; i >= 0; i--) {
aX[i] -= aY[i];
if (aX[i] < 0) {
if (i > 0) {
aX[i] += this.aBigInt;
aX[i - 1]--;
}
}
}
return aX;
};
getPi.prototype.isEmpty = function (aX) {
var empty = true;
for (var i = 0; i < aX.length; i++) {
if (aX[i]) {
empty = false;
break;
}
}
return empty;
};
getPi.prototype.getArcTan = function (iAng, n, aX) {
var iAng_squared = iAng * iAng;
var k = 3;
var sign = 0;
var arrAngle = [];
var aDivK = [];
aX = this.makeArr(n, aX, 0);
arrAngle = this.makeArr(n, arrAngle, 1);
arrAngle = this.Div(n, arrAngle, iAng, arrAngle);
aX = this.Add(n, aX, arrAngle);
while (!this.isEmpty(arrAngle)) {
arrAngle = this.Div(n, arrAngle, iAng_squared, arrAngle);
aDivK = this.Div(n, arrAngle, k, aDivK);
if (sign) {
aX = this.Add(n, aX, aDivK);
}
else {
aX = this.Sub(n, aX, aDivK);
}
k += 2;
sign = 1 - sign;
}
return aX;
};
getPi.prototype.makeArr = function (n1, arr, n2) {
for (var i = 0; i < n1; i++) {
arr[i] = null;
}
arr[0] = n2;
return arr;
};
return getPi;
}());
/*
tsc pi.ts && node pi.js
验证
3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679
计算所得(小数点后100位)
3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679
耗时: 0.001 秒
*/
var spaceSymbol = " ";
var app = new getPi(spaceSymbol);
var pi = app.getIt(100);
console.log(pi);把上述代码保存到文件 pi.js, 即可在 html 中引用，或者使用 Node.js 查看:
node pi.js
本文首发酷威普和51CTO,转载请注明。
本文出自 “13229289” 博客，转载请与作者联系！
以上是关于Machin formula /梅钦公式的主要内容，如果未能解决你的问题，请参考以下文章
python使用matplotlib可视化线图（line plot）可视化函数曲线并在标题title中显示函数的公式formula(add function formula into title)
罗德里格斯旋转公式（Rodrigues' rotation formula）推导