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,轉載請注明。