从圆面积计算π:-)

<input id="range" type="range" min="10" max="960" value="10" step="50" oninput="calcPi()"> <br> <div id="cont"></div> <script> function generateCircle(width) { var c = width/2; var delta = 1.0; var str = ""; var xCount = 0; for (var x=0; x <= width; x++) { for (var y = 0; y <= width; y++) { var d = Math.sqrt((x-c)*(x-c) + (y-c)*(y-c)); if (d > (width-1)/2) { str += '.'; } else { xCount++; str += 'o'; } str += "&nbsp;" } str += "\n"; } var pi = (xCount * 4) / (width * width); return [str, pi]; } function calcPi() { var e = document.getElementById("cont"); var width = document.getElementById("range").value; e.innerHTML = "<h4>Generating circle...</h4>"; setTimeout(function() { var circ = generateCircle(width); e.innerHTML = "<pre>" + "π = " + circ[1].toFixed(2) + "\n" + circ[0] +"</pre>"; }, 200); } calcPi(); </script>

下面是我在高中时学过的计算圆周率的技巧。

我之所以分享它，是因为我认为它足够简单，任何人都可以无限期地记住它，而且它教会了你“蒙特卡罗”方法的概念——这是一种统计方法，可以得到答案，这些答案不会立即通过随机过程演绎出来。

画一个正方形，在这个正方形内画一个象限(半圆的四分之一)(一个半径等于正方形边的象限，这样它就能尽可能多地填充正方形)

现在向正方形投掷飞镖，并记录飞镖落在何处——也就是说，在正方形内任意选择一个点。当然，它落在了正方形内部，但它落在半圆内部吗?记录这个事实。

重复此过程多次，你会发现半圆内的点数量与抛出的总数量之比为x。

由于正方形的面积是r乘以r，可以推导出半圆的面积是x乘以r乘以r(即x乘以r的平方)。因此x乘以4会得到。

这不是一个快速使用的方法。但这是蒙特卡罗方法的一个很好的例子。如果你环顾四周，你可能会发现许多超出你计算能力的问题都可以用这种方法来解决。

更好的方法

要获得标准常数(如pi)或标准概念的输出，我们应该首先使用所使用语言中可用的内置方法。它将以最快和最好的方式返回一个值。我正在使用python以最快的方式运行，以获得圆周率的值。

数学库的PI变量。数学库将变量pi存储为常数。

math_pi.py

import math
print math.pi

使用linux /usr/bin/time -v python math_pi.py的time工具运行脚本

输出:

Command being timed: "python math_pi.py"
User time (seconds): 0.01
System time (seconds): 0.01
Percent of CPU this job got: 91%
Elapsed (wall clock) time (h:mm:ss or m:ss): 0:00.03

用arccos的数学方法

acos_pi.py

import math
print math.acos(-1)

使用linux /usr/bin/time -v python acos_pi.py的time工具运行脚本

输出:

Command being timed: "python acos_pi.py"
User time (seconds): 0.02
System time (seconds): 0.01
Percent of CPU this job got: 94%
Elapsed (wall clock) time (h:mm:ss or m:ss): 0:00.03

使用BBP公式

bbp_pi.py

from decimal import Decimal, getcontext
getcontext().prec=100
print sum(1/Decimal(16)**k * 
          (Decimal(4)/(8*k+1) - 
           Decimal(2)/(8*k+4) - 
           Decimal(1)/(8*k+5) -
           Decimal(1)/(8*k+6)) for k in range(100))

使用linux /usr/bin/time -v python bbp_pi.py的time工具运行脚本

输出:

Command being timed: "python c.py"
User time (seconds): 0.05
System time (seconds): 0.01
Percent of CPU this job got: 98%
Elapsed (wall clock) time (h:mm:ss or m:ss): 0:00.06

因此，最好的方法是使用语言提供的内置方法，因为它们是获得输出的最快和最好的方法。在python中使用math.pi

实际上，有一本书专门介绍了π的快速计算方法:Jonathan和Peter Borwein写的《π和AGM》(在亚马逊上可以买到)。

我对年度股东大会和相关算法进行了相当多的研究:这非常有趣(尽管有时并非微不足道)。

请注意，要实现大多数现代算法来计算\pi，您将需要一个多精度算术库(GMP是一个很好的选择，尽管距离我上次使用它已经有一段时间了)。

最佳算法的时间复杂度为O(M(n)log(n))，其中M(n)是使用基于fft算法对两个n位整数(M(n)=O(n log(n) log(log(n)))相乘的时间复杂度，通常在计算\pi数字时需要fft算法，GMP中实现了该算法。

请注意，即使算法背后的数学可能并不简单，算法本身通常是几行伪代码，它们的实现通常非常简单(如果您选择不编写自己的多精度算术:-))。

从圆面积计算π:-)

<input id="range" type="range" min="10" max="960" value="10" step="50" oninput="calcPi()"> <br> <div id="cont"></div> <script> function generateCircle(width) { var c = width/2; var delta = 1.0; var str = ""; var xCount = 0; for (var x=0; x <= width; x++) { for (var y = 0; y <= width; y++) { var d = Math.sqrt((x-c)*(x-c) + (y-c)*(y-c)); if (d > (width-1)/2) { str += '.'; } else { xCount++; str += 'o'; } str += "&nbsp;" } str += "\n"; } var pi = (xCount * 4) / (width * width); return [str, pi]; } function calcPi() { var e = document.getElementById("cont"); var width = document.getElementById("range").value; e.innerHTML = "<h4>Generating circle...</h4>"; setTimeout(function() { var circ = generateCircle(width); e.innerHTML = "<pre>" + "π = " + circ[1].toFixed(2) + "\n" + circ[0] +"</pre>"; }, 200); } calcPi(); </script>

这是一个“经典”方法，非常容易实现。 这个在python(不是最快的语言)中的实现:

from math import pi
from time import time


precision = 10**6 # higher value -> higher precision
                  # lower  value -> higher speed

t = time()

calc = 0
for k in xrange(0, precision):
    calc += ((-1)**k) / (2*k+1.)
calc *= 4. # this is just a little optimization

t = time()-t

print "Calculated: %.40f" % calc
print "Constant pi: %.40f" % pi
print "Difference: %.40f" % abs(calc-pi)
print "Time elapsed: %s" % repr(t)

你可以在这里找到更多信息。

无论如何，在python中获得精确的圆周率值的最快方法是:

from gmpy import pi
print pi(3000) # the rule is the same as 
               # the precision on the previous code

下面是gpy pi方法的源代码，我认为在这种情况下，代码没有注释那么有用:

static char doc_pi[]="\
pi(n): returns pi with n bits of precision in an mpf object\n\
";

/* This function was originally from netlib, package bmp, by
 * Richard P. Brent. Paulo Cesar Pereira de Andrade converted
 * it to C and used it in his LISP interpreter.
 *
 * Original comments:
 * 
 *   sets mp pi = 3.14159... to the available precision.
 *   uses the gauss-legendre algorithm.
 *   this method requires time o(ln(t)m(t)), so it is slower
 *   than mppi if m(t) = o(t**2), but would be faster for
 *   large t if a faster multiplication algorithm were used
 *   (see comments in mpmul).
 *   for a description of the method, see - multiple-precision
 *   zero-finding and the complexity of elementary function
 *   evaluation (by r. p. brent), in analytic computational
 *   complexity (edited by j. f. traub), academic press, 1976, 151-176.
 *   rounding options not implemented, no guard digits used.
*/
static PyObject *
Pygmpy_pi(PyObject *self, PyObject *args)
{
    PympfObject *pi;
    int precision;
    mpf_t r_i2, r_i3, r_i4;
    mpf_t ix;

    ONE_ARG("pi", "i", &precision);
    if(!(pi = Pympf_new(precision))) {
        return NULL;
    }

    mpf_set_si(pi->f, 1);

    mpf_init(ix);
    mpf_set_ui(ix, 1);

    mpf_init2(r_i2, precision);

    mpf_init2(r_i3, precision);
    mpf_set_d(r_i3, 0.25);

    mpf_init2(r_i4, precision);
    mpf_set_d(r_i4, 0.5);
    mpf_sqrt(r_i4, r_i4);

    for (;;) {
        mpf_set(r_i2, pi->f);
        mpf_add(pi->f, pi->f, r_i4);
        mpf_div_ui(pi->f, pi->f, 2);
        mpf_mul(r_i4, r_i2, r_i4);
        mpf_sub(r_i2, pi->f, r_i2);
        mpf_mul(r_i2, r_i2, r_i2);
        mpf_mul(r_i2, r_i2, ix);
        mpf_sub(r_i3, r_i3, r_i2);
        mpf_sqrt(r_i4, r_i4);
        mpf_mul_ui(ix, ix, 2);
        /* Check for convergence */
        if (!(mpf_cmp_si(r_i2, 0) && 
              mpf_get_prec(r_i2) >= (unsigned)precision)) {
            mpf_mul(pi->f, pi->f, r_i4);
            mpf_div(pi->f, pi->f, r_i3);
            break;
        }
    }

    mpf_clear(ix);
    mpf_clear(r_i2);
    mpf_clear(r_i3);
    mpf_clear(r_i4);

    return (PyObject*)pi;
}

编辑:我在剪切和粘贴和缩进方面有一些问题，你可以在这里找到源代码。