梅钦公式( 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 实现:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 | /* * @Author: suifengtec * @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); |