什么是获得价值的最快方法?

任何语言都欢迎解决方案。 :-)我正在寻找最快的方式来获得?的价值,作为个人的挑战。更具体地说,我使用的方法不涉及使用 #define d常量,比如 M_PI ,或者对数字进行硬编码。

下面的程序测试了我所知道的各种方式。理论上,内联汇编版本是最快的选择,但显然不是可移植的;我已经将它作为基准来比较其他版本。在我的测试中,使用内置插件,GCC 4.2上的 4 * atan(1)版本是最快的,因为它将 atan(1)自动折叠为常量。指定 -fno-builtin 时, atan2(0,-1)版本是最快的。

这是主要的测试程序( pitimes.c ):

#include 
#include 
#include 

#define ITERS 10000000
#define TESTWITH(x) {                                                       \
    diff = 0.0;                                                             \
    time1 = clock();                                                        \
    for (i = 0; i < ITERS; ++i)                                             \
        diff += (x) - M_PI;                                                 \
    time2 = clock();                                                        \
    printf("%s\t=> %e, time => %f\n", #x, diff, diffclock(time2, time1));   \
}

static inline double
diffclock(clock_t time1, clock_t time0)
{
    return (double) (time1 - time0) / CLOCKS_PER_SEC;
}

int
main()
{
    int i;
    clock_t time1, time2;
    double diff;

    /* Warmup. The atan2 case catches GCC's atan folding (which would
     * optimise the ``4 * atan(1) - M_PI'' to a no-op), if -fno-builtin
     * is not used. */
    TESTWITH(4 * atan(1))
    TESTWITH(4 * atan2(1, 1))

#if defined(__GNUC__) && (defined(__i386__) || defined(__amd64__))
    extern double fldpi();
    TESTWITH(fldpi())
#endif

    /* Actual tests start here. */
    TESTWITH(atan2(0, -1))
    TESTWITH(acos(-1))
    TESTWITH(2 * asin(1))
    TESTWITH(4 * atan2(1, 1))
    TESTWITH(4 * atan(1))

    return 0;
}

和内联汇编的东西( fldpi.c ),指出它只适用于x86和x64系统:

double
fldpi()
{
    double pi;
    asm("fldpi" : "=t" (pi));
    return pi;
}

构建脚本构建我测试的所有配置( build.sh ):

#!/bin/sh
gcc -O3 -Wall -c           -m32 -o fldpi-32.o fldpi.c
gcc -O3 -Wall -c           -m64 -o fldpi-64.o fldpi.c

gcc -O3 -Wall -ffast-math  -m32 -o pitimes1-32 pitimes.c fldpi-32.o
gcc -O3 -Wall              -m32 -o pitimes2-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -fno-builtin -m32 -o pitimes3-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -ffast-math  -m64 -o pitimes1-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall              -m64 -o pitimes2-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall -fno-builtin -m64 -o pitimes3-64 pitimes.c fldpi-64.o -lm

除了在各种编译器标志之间进行测试之外(我还将32位与64位进行了比较,因为优化是不同的),我也尝试切换测试的顺序。尽管如此, atan2(0,-1)版本每次仍然是顶级的。

0
额外 编辑
意见: 17
@Zeus在这个特定的情况下,我的问题实际上是想成为一个“有趣的”微观优化问题(目前,它可能更适合编程拼图和代码高尔夫球),但”计算pi的最快方式“的一般前提似乎足以在这里保留这个问题。所以,在某个阶段,我可能会重新评估我是否应该接受最好的算法答案(可能是nlucaroni的答案),而不考虑它是否与微优化有关。
额外 作者 Chris Jester-Young,
@ HopelessN00b在我说的英语方言中,“优化”是拼写与“ s“,而不是”z“(发音为”zed“,顺便说一句,不是”zee“;-))。 (如果您查看评论历史,这不是我第一次也不得不恢复这种编辑。)
额外 作者 Chris Jester-Young,
nlucaroni的回答已经达到了100赞成(恭喜),所以它可能是一个很好的绿色勾号。请享用! (虽然,因为它是社区wiki和所有,它没有产生任何代表,所以甚至不知道如果nlucaroni会甚至注意到这一点。)
额外 作者 Chris Jester-Young,
额外 作者 Chris Jester-Young,
@erik:并不是所有的语言都有像 M_PI 这样的内置常量。我试图找到一种“权威”的方式来获得(理论上)跨多种语言(和/或其内置库)工作的pi(浮点)值。我目前首选的方法是使用 atan2(0,-1),但也许有更好的方法。
额外 作者 Chris Jester-Young,
@signonsridhar不,我们只讨论计算方法,它们在截断为双精度时给出与 M_PI 完全相同的结果。
额外 作者 Chris Jester-Young,
必须有一种方法可以在C ++元编程中实现。运行时间会非常好,但编译时间不会。
额外 作者 David Thornley,
只有一个解决方案比预先计算常量PI更快:预先计算公式中出现的所有值,例如,当需要周长时,可以预先计算2 * PI,而不是在运行时每次将PI乘以2。
额外 作者 ern0,
问题是:为什么你不想使用常量?例如或者由图书馆或自己定义?计算Pi是CPU周期的浪费,因为这个问题一遍又一遍地被解决了一些重要的数字,远远大于日常计算所需的数字
额外 作者 Tilo,
为什么你考虑使用atan(1)与使用M_PI不同?我会理解你为什么要这样做,如果你只使用算术运算,但是atan我没有看到这一点。
额外 作者 erikkallen,
理想情况下,这应该引起Mysticial的注意,因为他是计算pi的世界纪录保持者,其位数最多。 stackoverflow.com/questions/14283270/…
额外 作者 Team Upvote,
@ ChrisJester-Young不要。这不是我的意图。我最近开始更多地阅读规则......并且认为我会参与我访问的线索,提醒那些可能已经忘记了长时间框架的用户。我绝不想在这里当警察。道歉,如果我遇到粗鲁。
额外 作者 Zeus,
9801 /(1103?8)..给出了六位小数。这是计算PI的最快方法吗? = 3.14159273
额外 作者 signonsridhar,
@Chris Jester-Young我刚刚看到一个关于Ramanujan的视频,他给出了计算PI的方法。所以我只是分享它:>
额外 作者 signonsridhar,

20 答案

以下是我在高中学习计算pi的技术的一般描述。

我只是分享这一点,因为我认为这足够简单,任何人都可以无限期地记住它,再加上它会教给你“蒙特卡罗”方法的概念 - 这是统计方法得出的答案并不是立即出现通过随机过程可以推断出来。

绘制一个正方形,并在该正方形内部(半径等于正方形边的象限)内划一个象限(四分之一半圆),以便尽可能多地填充正方形)

现在在广场投掷一个飞镖,并记录它的落地位置 - 也就是说,在广场内的任何地方选择一个随机点。当然,它落在广场内,但它在半圆内吗?记录下这个事实。

重复这个过程很多次 - 你会发现半圈内的点数与抛出的总数的比值,称这个比率为x。

由于平方的面积是r乘以r,所以可以推断出半圆的面积是r乘以r的倍数(即x乘以r的平方)。因此x次4会给你pi。

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

0
额外
这是我们用来在学校的一个Java项目中计算Pi的方法。只是使用随机函数来提出x,y坐标,我们投掷的更多'飞镖'越接近我们来到的Pi。
额外 作者 Jeff Keslinke,

如果以最快的速度输入代码,请使用 golfscript 解决方案:

;''6666,-2%{2+.2/@*\/10.3??2*+}*`1000<~\;
0
额外

如上所述,蒙特卡罗方法应用了一些非常好的概念,但它显然不是最快的 - 不是一个远射,而是任何合理的用处。此外,这一切都取决于你正在寻找什么样的准确性。我知道的最快的pi是硬编码的数字。查看 PiPi [PDF] ,有很多公式。

这是一种快速收敛的方法(每次迭代约14位数)。当前最快的应用程序 PiFast 使用此公式与 FFT 。我只是写公式,因为代码很简单。这个公式几乎可以通过 Ramanujan和Chudnovsky发现找到。实际上他是如何计算数字的几十亿位数的 - 这不是一种无视的方法。由于我们正在分解阶乘因子,因此公式会很快溢出,因此延迟此类计算以删除项目将是有利的。

哪里,

以下是 Brent?Salamin算法。维基百科提到,当a和b足够接近时,(a + b)^ 2 / 4t将近似于pi。我不确定什么“足够接近”的意思,但从我的测试,一次迭代得到2digits,两个得到7,三个有15,当然这是双打,所以它可能有错误的基础上它的表示和'真实“的计算可能更准确。

let pi_2 iters =
    let rec loop_ a b t p i =
        if i = 0 then a,b,t,p
        else
            let a_n = (a +. b) /. 2.0 
            and b_n = sqrt (a*.b)
            and p_n = 2.0 *. p in
            let t_n = t -. (p *. (a -. a_n) *. (a -. a_n)) in
            loop_ a_n b_n t_n p_n (i - 1)
    in 
    let a,b,t,p = loop_ (1.0) (1.0 /. (sqrt 2.0)) (1.0/.4.0) (1.0) iters in
    (a +. b) *. (a +. b) /. (4.0 *. t)

最后,一些pi高尔夫(800位数字)怎么样? 160个字符!

int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}
0
额外
假设你试图自己实现第一个,sqr(k3)不会成为问题吗?我很确定这会导致一个不合理的数字,你将不得不估计(IIRC,所有不是整数的根都是非理性的)。如果您使用无限精度算术,但其平方根是一个交易断路器,其他所有东西都很直观。第二个包括一个sqrt。
额外 作者 Bill K,
根据我的经验,“足够接近”通常意味着涉及泰勒级数逼近。
额外 作者 Stephen,

我真的很喜欢这个程序,它通过查看它自己的区域近似于pi :-)

IOCCC 1988: westley.c

#define _ -F<00||--F-OO--;
int F=00,OO=00;main(){F_OO();printf("%1.3f\n",4.*-F/OO/OO);}F_OO()
{
            _-_-_-_
       _-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
        _-_-_-_-_-_-_-_
            _-_-_-_
}
0
额外
如果用_F <00 || --F-OO替换_,应该更容易遵循:-)
额外 作者 Pat,
这个程序在1998年非常棒,但因为现代预处理程序比较自由,在宏扩展周围插入空格以防止这样的事情发挥作用,所以这个程序被破坏了。不幸的是,这是一件遗物。
额外 作者 Chris Lutz,
它在这里打印0.25 -
额外 作者 Johannes Schaub - litb,
@Pat如果你为什么编辑它,那是因为我在LQP队列中看到了这个答案 stackoverflow.com/review/low -quality-posts / 16750528 ,因此为了避免删除,我在链接中添加了代码。
额外 作者 Petter Friberg,
- traditional-cpp 传递给 cpp 以获得预期的行为。
额外 作者 Nietzche-jou,
或者,如果您将_替换为“if(前一个字符是' - '){OO--;} F--;”
额外 作者 FryGuy,

早在过去,由于字体小,缓存或不存在的浮点操作,我们曾经这样做过:

/* Return approximation of n * PI; n is integer */
#define pi_times(n) (((n) * 22) / 7)

对于不需要很高精度的应用(例如视频游戏),这非常快速并且足够准确。

0
额外
为了更准确地使用 355/113 。非常准确的数字大小涉及。
额外 作者 David Thornley,
出于好奇: 22/73 + 1/7
额外 作者 Agnius Vasiliauskas,

您可以使用 BBP公式计算第n位数 - 以2为基数(或16 ) - 无需先打扰先前的n-1数字:)

0
额外

这个版本(在Delphi中)没有特别之处,但它至少比版本更快Nick Hodge在他的博客上发布了 :)。在我的机器上,执行十亿次迭代大约需要16秒,给出的值为 3.14159265 25879(准确的部分用粗体表示)。

program calcpi;

{$APPTYPE CONSOLE}

uses
  SysUtils;

var
  start, finish: TDateTime;

function CalculatePi(iterations: integer): double;
var
  numerator, denominator, i: integer;
  sum: double;
begin
  {
  PI may be approximated with this formula:
  4 * (1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 .......)
  //}
  numerator := 1;
  denominator := 1;
  sum := 0;
  for i := 1 to iterations do begin
    sum := sum + (numerator/denominator);
    denominator := denominator + 2;
    numerator := -numerator;
  end;
  Result := 4 * sum;
end;

begin
  try
    start := Now;
    WriteLn(FloatToStr(CalculatePi(StrToInt(ParamStr(1)))));
    finish := Now;
    WriteLn('Seconds:' + FormatDateTime('hh:mm:ss.zz',finish-start));
  except
    on E:Exception do
      Writeln(E.Classname, ': ', E.Message);
  end;
end.
0
额外

如果这篇文章为真,那么 Bellard 创建的算法可能是最快速的可用算法之一。他使用台式电脑将pi创建为2.7万亿位数!

...and he has published his work here

做好工作Bellard,你是先驱!

http://www.theregister.co.uk/2010/01/06/ very_long_pi /

0
额外
Bellard的开创方式有很多种,首先是LZEXE,很可能是第一个可执行压缩器(想想UPX的作用,然后回溯到80年代),当然,QEMU和FFMPEG都被广泛使用。哦,他的IOCCC入口.... :-P
额外 作者 Chris Jester-Young,

如果您愿意使用近似值, 355/113 适用于6位十进制数字,并且具有可用于整数表达式的附加优点。现在这并不重要,因为“浮点数学协处理器”不再具有任何意义,但它一度相当重要。

0
额外

为了完整起见,一个C ++模板版本,对于一个优化版本,它将在编译时计算PI,并将内联到一个单一值。

#include 

template
struct sign
{
    enum {value = (I % 2) == 0 ? 1 : -1};
};

template
struct pi_calc
{
    inline static double value ()
    {
        return (pi_calc::value() + pi_calc::value ()) / 2.0;
    }
};

template
struct pi_calc<0, J>
{
    inline static double value ()
    {
        return (sign::value * 4.0) / (2.0 * J + 1.0) + pi_calc<0, J-1>::value ();
    }
};


template<>
struct pi_calc<0, 0>
{
    inline static double value ()
    {
        return 4.0;
    }
};

template
struct pi
{
    inline static double value ()
    {
        return pi_calc::value ();
    }
};

int main ()
{
    std::cout.precision (12);

    const double pi_value = pi<10>::value ();

    std::cout << "pi ~ " << pi_value << std::endl;

    return 0;
}

Note for I > 10, optimised builds can be slow, likewise for non-optimised runs. For 12 iterations I believe there are around 80k calls to value() (in the absence of memoisation).

0
额外
@sebastião-miranda 莱布尼茨公式,平均加速度提高了收敛性。 pi_calc <0,J> 计算公式中的每个连续项,非特定的 pi_calc 计算平均值。
额外 作者 jon-hanson,
那么,这是准确的9dp。你是在反对某件事还是只是在做一个观察?
额外 作者 jon-hanson,
我运行这个并获得“pi〜3.14159265383”
额外 作者 maxwellb,
这里用来计算PI的算法的名称是什么?
额外 作者 Sebastião Miranda,

计算?来自圈子区域:-)

<div class="snippet" data-lang="js" data-hide="false" data-console="true" data-babel="false"> <div class="snippet-code">

<input id="range" type="range" min="10" max="960" value="10" step="50" oninput="calcPi()">

<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 += " " } 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 = "

Generating circle...

";
    setTimeout(function() {
        var circ = generateCircle(width);
        e.innerHTML  = "
" + "? = " + circ[1].toFixed(2) + "\n" + circ[0] +"
";
    }, 200);
}
calcPi();
</script>
</div> </div>

0
额外

以下答案正是如何以最快的方式实现这一点 - 以最少的计算工作量。即使你不喜欢答案,你也必须承认这确实是获得PI价值的最快方法。

获取Pi值的最快方式是:

  1. 选择了您最喜欢的编程语言
  2. 加载它的数学库
  3. ,并发现Pi已经在那里定义了!准备好使用它..

如果你手边没有数学库。

SECOND FASTEST 方式(更通用的解决方案)是:

在互联网上查找Pi,例如这里:

http://www.eveandersson.com/pi/digits/1000000 (1 million digits .. what's your floating point precision? )

或在这里:

http://3.141592653589793238462643383279502884197169399375105820974944592.com/

或在这里:

http://en.wikipedia.org/wiki/Pi

无论您想要使用什么精度算法,您都可以快速找到需要的数字,通过定义一个常数,可以确保您不会浪费宝贵的CPU时间。

这不仅仅是一个幽默的回答,实际上,如果有人愿意在实际应用中计算出Pi的价值,那将会是一个相当大的CPU时间浪费,不是吗?至少我没有看到一个真正的应用程序试图重新计算这个。

尊敬的主持人:请注意,OP询问:最快的方式以获得PI的价值“

0
额外

如果你想计算值的近似值, (出于某种原因),你应该尝试一个二进制提取算法。 Bellard's 改进 BBP 给出O(N ^ 2)中的PI。


如果你想获得值的近似值,做计算,然后:

PI = 3.141592654

当然,这只是一个近似值,并不完全准确。它超过了0.00000000004102。 (四十万亿,大约四亿分之一)。


如果你想用<�?>做数学,那么给自己一个铅笔和纸或者一个计算机代数包,并且使用?的确切值。

如果你真的想要一个公式,这个很有趣:

? = - i ln(-1)

0
额外
你的公式取决于你如何在复杂平面中定义ln。它必须在复平面中沿着一条线不连续,并且该线很常见为负实轴。
额外 作者 erikkallen,

刚刚遇到这个应该在这里完整的一个:

计算PI中的PI

它具有相当不错的属性,可以提高精度,使程序更大。

Here's some insight into the language itself

0
额外

用D编译时计算PI。

(复制自 DSource.org

/** Calculate pi at compile time
 *
 * Compile with dmd -c pi.d
 */
module calcpi;

import meta.math;
import meta.conv;

/** real evaluateSeries!(real x, real metafunction!(real y, int n) term)
 *
 * Evaluate a power series at compile time.
 *
 * Given a metafunction of the form
 *  real term!(real y, int n),
 * which gives the nth term of a convergent series at the point y
 * (where the first term is n==1), and a real number x,
 * this metafunction calculates the infinite sum at the point x
 * by adding terms until the sum doesn't change any more.
 */
template evaluateSeries(real x, alias term, int n=1, real sumsofar=0.0)
{
  static if (n>1 && sumsofar == sumsofar + term!(x, n+1)) {
     const real evaluateSeries = sumsofar;
  } else {
     const real evaluateSeries = evaluateSeries!(x, term, n+1, sumsofar + term!(x, n));
  }
}

/*** Calculate atan(x) at compile time.
 *
 * Uses the Maclaurin formula
 *  atan(z) = z - z^3/3 + Z^5/5 - Z^7/7 + ...
 */
template atan(real z)
{
    const real atan = evaluateSeries!(z, atanTerm);
}

template atanTerm(real x, int n)
{
    const real atanTerm =  (n & 1 ? 1 : -1) * pow!(x, 2*n-1)/(2*n-1);
}

/// Machin's formula for pi
/// pi/4 = 4 atan(1/5) - atan(1/239).
pragma(msg, "PI = " ~ fcvt!(4.0 * (4*atan!(1/5.0) - atan!(1/239.0))) );
0
额外
不幸的是,切线是以pi为基础的,这种计算有些无效。
额外 作者 Grant Johnson,

实际上,有一本书是专门用于计算\ pi:'Pi和AGM'的 fast 方法的一本书,由Jonathan和Peter Borwein( rel =”noreferrer“title =”在Amazon上可用“)。

我研究了AGM和相关算法:它非常有趣(尽管有时并不重要)。

请注意,要实现大多数现代算法来计算\ pi,您需要一个多精度算术库(

GMP 是一个不错的选择,尽管自从我上次使用它以来已经有一段时间了)。

最佳算法的时间复杂度为O(M(n)log(n)),其中M(n)是两个n位整数相乘的时间复杂度(M(n)= O(n log(n)log(log(n)))使用基于FFT的算法,这些算法通常在计算\ pi的数字时需要,并且这种算法在GMP中实现)。

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

0
额外

这是一种“经典”方法,非常容易实现。 这个实现,在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 "Costant pi: %.40f" % pi
print "Difference: %.40f" % abs(calc-pi)
print "Time elapsed: %s" % repr(t)

您可以在此处找到更多信息。

无论如何,在python中获得一个精确的,尽可能多的pi值的最快方法是:

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

这里是gmpy 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;
}

EDIT: I had some problem with cut and paste and identation, anyway you can find the source here.

0
额外

克里斯上面布兰特的方法非常好;布伦特通常是任意精度算术领域的巨人。

如果你想要的是第N位,那就是着名的 BBP公式 在十六进制中很有用

0
额外
布伦特方法不是由我发布的;它由Andrea发布,我恰好是编辑帖子的最后一个人。 :-)但我同意,该帖子值得赞赏。
额外 作者 Chris Jester-Young,

双打:

4.0 * (4.0 * Math.Atan(0.2) - Math.Atan(1.0 / 239.0))

这将精确到小数点后14位,足以填充双精度(不准确可能是因为圆弧切线中的其余小数点被截断)。

还有Seth,它是3.14159265358979323846 3 ,而不是64。

0
额外

使用类似Machin的公式

176 * arctan (1/57) + 28 * arctan (1/239) - 48 * arctan (1/682) + 96 * arctan(1/12943) 

[; \left( 176 \arctan \frac{1}{57} + 28 \arctan \frac{1}{239} - 48 \arctan \frac{1}{682} + 96 \arctan \frac{1}{12943}\right) ;], for you TeX the World people.

在Scheme中实现,例如:

(*(atan(/ 1 682))))*(* 96(atan(1/157))) (atan(/ 1 12943))))

0
额外