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)公式，但梅钦公式至今仍然是计算 π 值的主要公式。

梅钦公式是这样的：

Machin formula

或者写为：

Machin formula

梅钦公式是格里高利/莱布尼茨计算公式的变体，但是更实用，它的收敛速度显著增加，这使得它成为了更实用的计算π的方法。

实现

JavaScript 实现:

/*
* @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);

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);

π 的简史

简介

实现

相关文章