Nomar记录一段历史
获得 π 值的最快方法是什么?

我正在寻找获得 π 值的最快方法,作为个人挑战。更具体地说,我使用的方法不涉及使用#define常量(如M_PI )或对数字进行硬编码。

下面的程序测试了我所知道的各种方法。理论上,内联汇编版本是最快的选择,但显然不可移植。我已将其作为基准与其他版本进行比较。在我的测试中,使用内置函数, 4 * atan(1)版本在 GCC 4.2 上速度最快,因为它会自动将atan(1)折叠成一个常量。指定-fno-builtinatan2(0, -1)版本最快。

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

 #include <math.h> #include <stdio.h> #include <time.h> #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; }

以及仅适用于 x86 和 x64 系统的内联汇编( fldpi.c ):

 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)版本仍然每次都名列前茅。

刚刚遇到了这个应该在这里的完整性:在 Piet 中计算 PI它具有相当好的属性,可以提高精度使程序更大。是对语言本身的一些见解
为什么您认为使用 atan(1) 与使用 M_PI 不同?如果您只使用算术运算,我会理解为什么要这样做,但是对于 atan 我不明白这一点。
@erik:并非所有语言都有像M_PI这样的内置常量。我试图找到一种“权威”方法来获得 pi 的(浮点)值,该值(理论上)适用于各种语言(和/或其内置库)。我目前首选的方法是使用atan2(0, -1) ,但也许有更好的方法。
Bellard 在许多方面开创了先河……首先是 LZEXE,很可能是第一个可执行压缩器(想想 UPX 的作用,然后回到 80 年代),当然现在,QEMU 和 FFMPEG 都被广泛使用.哦,还有他的 IOCCC 条目....:-P
问题是:为什么你希望使用一个常数?例如,要么由图书馆定义,要么由您自己定义?计算 Pi 是对 CPU 周期的浪费,因为这个问题已经被一遍又一遍地解决,其有效位数远大于日常计算所需的位数
只有一种解决方案比预先计算常数 PI 更快:预先计算公式中出现的所有值,例如,当需要周长时,您可以预先计算 2*PI,而不是在运行时每次将 PI 乘以 2。
@ChrisJester-Young 不。这不是意图.. 我最近开始阅读更多规则.. 并认为我会在我访问的线程中尽自己的一份力量来提醒可能已经忘记了很长一段时间的用户。我绝对不想在这里当警察。如果我遇到粗鲁,我深表歉意。
@Zeus 在这种特定情况下,我的问题实际上是一个“有趣”的微优化问题(现在,这可能更适合编程拼图和代码高尔夫),但是“最快方式”的一般前提计算圆周率”似乎足够有用,可以将这个问题留在这里。因此,在某个阶段,我可能会重新评估我是否应该只接受最佳算法答案(可能是 nlucaroni 的答案),而不考虑它是否与微优化有关。
@HopelessN00b 在我说的英语方言中,“优化”拼写为“s”,而不是“z”(发音为“zed”,顺便说一句,而不是“zee”;-))。 (如果您查看评论历史记录,这也不是我第一次不得不恢复此类编辑。)
nlucaroni 的回答已经达到了 100 票(恭喜),所以它可能是一个很好的点。享受! (不过,由于它是社区 wiki 和所有内容,它不会产生任何代表,所以甚至不确定 nlucaroni 是否会注意到这一点。)
9801/(1103√8)..给出六位小数..这是计算PI的最快方法吗? = 3.14159273
@signonsridhar 不,我们只讨论在截断为双精度时给出与M_PI完全相同结果的计算方法。
@Chris Jester-Young 好吧,我刚刚看到了一个关于 Ramanujan 的视频,他用这种方法计算 PI。所以我只是分享了它:>
您可以从学习如何获得 π 的公式开始,然后尝试将其写在 C 文件中。使用Brent-Salamin 公式,您可以尝试找出其算法并将其转换为程序。我还不知道该怎么做,因为我不明白一些事情......

22个回答

这是我在高中学到的计算圆周率的技术的一般描述。

我只是分享这个,因为我认为它很简单,任何人都可以无限期地记住它,而且它教会了你“蒙特卡罗”方法的概念——这是一种统计方法,可以得出不会立即出现的答案可以通过随机过程推导出来。

画一个正方形,并在该正方形内刻出一个象限(半圆的四分之一)(半径等于正方形边长的象限,因此它尽可能多地填充正方形)

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

多次重复这个过程——你会发现半圆内的点数与抛出的总数之比,​​称这个比率为x。

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

这不是一个快速的使用方法。但它是蒙特卡罗方法的一个很好的例子。如果您环顾四周,您可能会发现许多超出您计算技能的问题都可以通过此类方法解决。

这是我们在学校的一个java项目中用来计算Pi的方法。只是使用随机器来计算 x,y 坐标,我们投掷的“飞镖”越多,距离 Pi 越近。

如前所述,蒙特卡罗方法应用了一些伟大的概念,但显然,它不是最快的,也不是远景,也不是任何合理的衡量标准。此外,这完全取决于您正在寻找什么样的准确性。我所知道的最快的 π 是数字硬编码的那个。看PiPi[PDF] ,有很多公式。

这是一种快速收敛的方法——每次迭代大约 14 位数字。 PiFast是当前最快的应用程序,使用此公式和 FFT。我只会写公式,因为代码很简单。这个公式几乎被拉马努金发现,被楚德诺夫斯基发现。这实际上是他计算出的数十亿位数字的方式——所以这不是一种可以无视的方法。该公式将很快溢出,并且由于我们正在除法阶乘,因此延迟此类计算以删除项将是有利的。

在此处输入图片说明

在此处输入图片说明

在哪里,

在此处输入图片说明

下面是Brent-Salamin 算法。维基百科提到,当ab “足够接近”时, (a + b)² / 4t将是 π 的近似值。我不确定“足够接近”是什么意思,但从我的测试来看,一次迭代得到 2 位数,两次得到 7,三次得到 15,当然这是双打,所以根据它的表示它可能有错误和真正的计算可能更准确。

 let pi_2 iters = let rec loop_ abtpi = 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)

最后,来点圆周率高尔夫(800 位数)怎么样? 160个字!

 int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;bc;)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);}
假设您正在尝试自己实现第一个,那么 sqr(k3) 不会有问题吗?我很确定它最终会得到一个你必须估计的无理数(IIRC,所有不是整数的根都是无理数)。如果您使用无限精度算术,其他一切看起来都非常简单,但平方根是一个交易破坏者。第二个也包括一个 sqrt。
根据我的经验,“足够接近”通常意味着涉及泰勒级数近似。

如果最快的意思是输入代码最快,这里是Golfscript解决方案:

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

实际上有一整本书专门(除其他外)用于快速计算 \pi 的方法:'Pi and the AGM',作者 Jonathan 和 Peter Borwein( 可在亚马逊上获得)。

我研究了很多 AGM 和相关算法:它很有趣(虽然有时很重要)。

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

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

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

BBP 公式允许您计算第 n 个数字 - 以 2(或 16)为基数 - 甚至不必先考虑前面的 n-1 个数字:)

我真的很喜欢这个程序,因为它通过查看自己的面积来近似 π。

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() { _-_-_-_ _-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_ _-_-_-_ }
如果您将 _ 替换为 -F<00||--F-OO-- 应该更容易遵循:-)
或者,如果您将 _ 替换为“if(前一个字符是 '-'){ OO--; } F--;”
这个程序在 1998 年很棒,但被破坏了,因为现代预处理器更加自由地在宏扩展周围插入空格以防止这样的事情工作。不幸的是,这是一个遗物。
--traditional-cpp传递给cpp以获得预期的行为。
@Pat 如果你更清楚为什么我编辑它是因为我在 LQP 队列stackoverflow.com/review/low-quality-posts/16750528 中看到了这个答案,因此为了避免删除,我在答案的链接中添加了代码。

在编译时用 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))) );
不幸的是,正切是基于 pi 的反正切,这在某种程度上使此计算无效。

这是一个“经典”的方法,很容易实现。 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中获得精确的你想要的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 rp brent), in analytic computational * complexity (edited by jf 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; }

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

这个版本(在 Delphi 中)没什么特别的,但它至少比Nick Hodge 在他的博客上发布的版本要快 :)。在我的机器上,执行 10 亿次迭代大约需要 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.

回到过去,由于字长小,浮点运算缓慢或不存在,我们过去常常做这样的事情:

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

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

为了更准确,请使用355 / 113 。对于所涉及的数字的大小非常准确。

圆周率正好是 3! [教授。弗林克(辛普森一家)]

笑话,但这是 C# 中的一个(需要 .NET-Framework)。

 using System; using System.Text; class Program { static void Main(string[] args) { int Digits = 100; BigNumber x = new BigNumber(Digits); BigNumber y = new BigNumber(Digits); x.ArcTan(16, 5); y.ArcTan(4, 239); x.Subtract(y); string pi = x.ToString(); Console.WriteLine(pi); } } public class BigNumber { private UInt32[] number; private int size; private int maxDigits; public BigNumber(int maxDigits) { this.maxDigits = maxDigits; this.size = (int)Math.Ceiling((float)maxDigits * 0.104) + 2; number = new UInt32[size]; } public BigNumber(int maxDigits, UInt32 intPart) : this(maxDigits) { number[0] = intPart; for (int i = 1; i < size; i++) { number[i] = 0; } } private void VerifySameSize(BigNumber value) { if (Object.ReferenceEquals(this, value)) throw new Exception("BigNumbers cannot operate on themselves"); if (value.size != this.size) throw new Exception("BigNumbers must have the same size"); } public void Add(BigNumber value) { VerifySameSize(value); int index = size - 1; while (index >= 0 && value.number[index] == 0) index--; UInt32 carry = 0; while (index >= 0) { UInt64 result = (UInt64)number[index] + value.number[index] + carry; number[index] = (UInt32)result; if (result >= 0x100000000U) carry = 1; else carry = 0; index--; } } public void Subtract(BigNumber value) { VerifySameSize(value); int index = size - 1; while (index >= 0 && value.number[index] == 0) index--; UInt32 borrow = 0; while (index >= 0) { UInt64 result = 0x100000000U + (UInt64)number[index] - value.number[index] - borrow; number[index] = (UInt32)result; if (result >= 0x100000000U) borrow = 0; else borrow = 1; index--; } } public void Multiply(UInt32 value) { int index = size - 1; while (index >= 0 && number[index] == 0) index--; UInt32 carry = 0; while (index >= 0) { UInt64 result = (UInt64)number[index] * value + carry; number[index] = (UInt32)result; carry = (UInt32)(result >> 32); index--; } } public void Divide(UInt32 value) { int index = 0; while (index < size && number[index] == 0) index++; UInt32 carry = 0; while (index < size) { UInt64 result = number[index] + ((UInt64)carry << 32); number[index] = (UInt32)(result / (UInt64)value); carry = (UInt32)(result % (UInt64)value); index++; } } public void Assign(BigNumber value) { VerifySameSize(value); for (int i = 0; i < size; i++) { number[i] = value.number[i]; } } public override string ToString() { BigNumber temp = new BigNumber(maxDigits); temp.Assign(this); StringBuilder sb = new StringBuilder(); sb.Append(temp.number[0]); sb.Append(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.CurrencyDecimalSeparator); int digitCount = 0; while (digitCount < maxDigits) { temp.number[0] = 0; temp.Multiply(100000); sb.AppendFormat("{0:D5}", temp.number[0]); digitCount += 5; } return sb.ToString(); } public bool IsZero() { foreach (UInt32 item in number) { if (item != 0) return false; } return true; } public void ArcTan(UInt32 multiplicand, UInt32 reciprocal) { BigNumber X = new BigNumber(maxDigits, multiplicand); X.Divide(reciprocal); reciprocal *= reciprocal; this.Assign(X); BigNumber term = new BigNumber(maxDigits); UInt32 divisor = 1; bool subtractTerm = true; while (true) { X.Divide(reciprocal); term.Assign(X); divisor += 2; term.Divide(divisor); if (term.IsZero()) break; if (subtractTerm) this.Subtract(term); else this.Add(term); subtractTerm = !subtractTerm; } } }

我总是使用acos(-1) ,而不是将 pi 定义为常量。

cos(-1) 还是 acos(-1)? :-P 那(后者)是我原始代码中的测试用例之一。它是我的首选(以及 atan2(0, -1),它确实与 acos(-1) 相同,只是 acos 通常是根据 atan2 实现的),但有些编译器针对 4 * atan(1) 进行了优化!

如果您愿意使用近似值, 355 / 113适用于 6 位十进制数字,并且具有可用于整数表达式的额外优势。如今,这并不重要,因为“浮点数学协处理器”不再具有任何意义,但它曾经非常重要。

出于完整性考虑,C++ 模板版本,对于优化构建,将在编译时计算 PI 的近似值,并将内联到单个值。

 #include <iostream> template<int I> struct sign { enum {value = (I % 2) == 0 ? 1 : -1}; }; template<int I, int J> struct pi_calc { inline static double value () { return (pi_calc<I-1, J>::value () + pi_calc<I-1, J+1>::value ()) / 2.0; } }; template<int J> struct pi_calc<0, J> { inline static double value () { return (sign<J>::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<int I> struct pi { inline static double value () { return pi_calc<I, I>::value (); } }; int main () { std::cout.precision (12); const double pi_value = pi<10>::value (); std::cout << "pi ~ " << pi_value << std::endl; return 0; }

请注意,对于 I > 10,优化构建可能会很慢,对于非优化运行也是如此。对于 12 次迭代,我相信大约有 80k 次调用 value()(在没有记忆的情况下)。

我运行它并得到“pi ~ 3.14159265383”
嗯,这对 9dp 来说是准确的。你是在反对某事还是只是在做一个观察?
这里用来计算 PI 的算法的名称是什么?
@sebastião-miranda Leibniz 的公式,平均加速度可提高收敛性。 pi_calc<0, J>计算公式中的每个连续项,非专业的pi_calc<I, J>计算平均值。

如果您想计算π 值的近似值(出于某种原因),您应该尝试二进制提取算法。 BellardBBP改进给出了 O(N^2) 中的 PI。


如果要获得π 值的近似值来进行计算,则:

 PI = 3.141592654

当然,这只是一个近似值,并不完全准确。它比 0.00000000004102 多一点。 (四万亿分之四,大约4 / 10,000,000,000 )。


如果你想用 π 做数学,那么给自己准备一支纸笔或一个计算机代数包,并使用 π 的精确值 π。

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

π = - i ln(-1)

您的公式取决于您如何在复平面中定义 ln。它必须沿着复平面中的一条线不连续,并且这条线是负实轴是很常见的。

双打:

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

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

同样是赛斯,它是 3.14159265358979323846 3 ,而不是 64。

使用类似机器的公式

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 中实现,例如:

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

下面准确地回答了如何以最快的方式执行此操作 - 使用最少的计算工作。即使你不喜欢这个答案,你也不得不承认,这确实是获得 PI 价值的最快方式。

获取 Pi 值的最快方法是:

1) 选择你最喜欢的编程语言 2) 加载它的数学库 3) 发现 Pi 已经在那里定义了——可以使用了!

如果您手头没有数学库..

第二种最快的方式(更通用的解决方案)是:

在 Internet 上查找 Pi,例如:

http://www.eveandersson.com/pi/digits/1000000(100万位数字......你的浮点精度是多少?)

或在这里:

http://3.141592653589793238462643383279502884197169399375105820974944592.com/

或在这里:

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

为您想要使用的任何精度算术找到所需的数字真的很快,并且通过定义一个常量,您可以确保不会浪费宝贵的 CPU 时间。

这不仅是一个部分幽默的答案,而且在现实中,如果有人继续计算实际应用程序中 Pi 的值......那将是相当大的 CPU 时间浪费,不是吗?至少我没有看到尝试重新计算这个的真正应用程序。

尊敬的版主:请注意OP问:“获取PI值的最快方式”

亲爱的 Tilo:请注意 OP 说:“我正在寻找获得 π 值的最快方法,作为个人挑战。更具体地说,我正在使用不涉及使用 #define 常量(如 M_PI)的方法,或硬编码中的数字
亲爱的@Max:请注意,OP我回答编辑了他们的原始问题 - 这几乎不是我的错;) 我的解决方案仍然是最快的方法,并且以任何所需的浮点精度和没有 CPU 周期优雅地解决了问题:)
哦对不起,我没有意识到。只是想一想,硬编码常量的精度不会比计算 pi 的精度低吗?我想这取决于它是什么语言以及创建者愿意将所有数字放入 :-)

从圆面积计算 π :-)

 <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((xc)*(xc) + (yc)*(yc)); 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>

更好的方法

要获得pi或标准概念等标准常量的输出,我们应该首先使用您所使用的语言中可用的内置方法。它将以最快和最好的方式返回一个值。我正在使用python运行获取pi值的最快方法。

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

math_pi.py

 import math print math.pi

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

输出:

 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
  • 使用arc cos数学方法

acos_pi.py

 import math print math.acos(-1)

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

输出:

 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_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时间实用程序运行脚本

输出:

 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

如果您不介意执行平方根和几个逆运算,那么 Chudnovsky 算法非常快。它仅在 2 次迭代中收敛到双精度。

 /* Chudnovsky algorithm for computing PI */ #include <iostream> #include <cmath> using namespace std; double calc_PI(int K=2) { static const int A = 545140134; static const int B = 13591409; static const int D = 640320; const double ID3 = 1./ (double(D)*double(D)*double(D)); double sum = 0.; double b = sqrt(ID3); long long int p = 1; long long int a = B; sum += double(p) * double(a)* b; // 2 iterations enough for double convergence for (int k=1; k<K; ++k) { // A*k + B a += A; // update denominator b *= ID3; // p = (-1)^k 6k! / 3k! k!^3 p *= (6*k)*(6*k-1)*(6*k-2)*(6*k-3)*(6*k-4)*(6*k-5); p /= (3*k)*(3*k-1)*(3*k-2) * k*k*k; p = -p; sum += double(p) * double(a)* b; } return 1./(12*sum); } int main() { cout.precision(16); cout.setf(ios::fixed); for (int k=1; k<=5; ++k) cout << "k = " << k << " PI = " << calc_PI(k) << endl; return 0; }

结果:

 k = 1 PI = 3.1415926535897341 k = 2 PI = 3.1415926535897931 k = 3 PI = 3.1415926535897931 k = 4 PI = 3.1415926535897931 k = 5 PI = 3.1415926535897931

基本上是回形针优化器答案的 C 版本,并且更加简化:

 #include <stdio.h> #include <math.h> double calc_PI(int K) { static const int A = 545140134; static const int B = 13591409; static const int D = 640320; const double ID3 = 1.0 / ((double) D * (double) D * (double) D); double sum = 0.0; double b = sqrt(ID3); long long int p = 1; long long int a = B; sum += (double) p * (double) a * b; for (int k = 1; k < K; ++k) { a += A; b *= ID3; p *= (6 * k) * (6 * k - 1) * (6 * k - 2) * (6 * k - 3) * (6 * k - 4) * (6 * k - 5); p /= (3 * k) * (3 * k - 1) * (3 * k - 2) * k * k * k; p = -p; sum += (double) p * (double) a * b; } return 1.0 / (12 * sum); } int main() { for (int k = 1; k <= 5; ++k) { printf("k = %i, PI = %.16f\n", k, calc_PI(k)); } }

但是为了更简单,这个算法采用了Chudnovsky的公式,如果你不是真正理解代码,我可以完全简化。

总结:我们将获得一个从 1 到 5 的数字并将其添加到我们将用来获得 PI 的函数中。然后给你3个号码:545140134(A)、13591409(B)、640320(D)。然后我们将使用 D 作为一个double乘以 3 倍到另一个double (ID3)。然后我们将 ID3 的平方根带入另一个double精度数 (b) 并分配 2 个数字:1 (p),B (a) 的值。请注意,C 不区分大小写。然后将通过乘以 p、a 和 b 的值来创建double (总和),所有这些都在double精度中。然后循环直到给函数给出的数字将开始并将 A 的值与 a 相加,b 的值乘以 ID3,p 的值将乘以多个值,我希望您能理解,并且还除以多个值作为好。总和将再次与 p、a 和 b 相加,循环将重复,直到循环数的值大于或等于 5。随后,总和乘以 12 并由函数返回,结果为PI。

好吧,那很长,但我想你会明白的……

随机文章