Какой самый быстрый способ получить значение π?

Мне очень нравится эта программа, потому что она приближается к π, глядя на ее собственную область.

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()
{
            _-_-_-_
       _-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
        _-_-_-_-_-_-_-_
            _-_-_-_
}

Вот общее описание метода вычисления pi, который я изучил в старшей школе.

Я разделяю это только потому, что считаю, что это достаточно просто, что каждый может запомнить его бесконечно, плюс он учит вас концепции методов "Монте-Карло", которые являются статистическими методами получения ответов, которые не сразу кажутся выводимыми через случайные процессы.

Нарисуйте квадрат и вставьте квадрат (четверть полукруга) внутри этого квадрата (квадрант с радиусом, равным стороне квадрата, поэтому он заполняет как можно больше квадрата)

Теперь бросьте дротик на квадрат и записывайте, где он приземляется, то есть выберите случайную точку в любом месте квадрата. Конечно, он приземлился на площади, но находится ли он в полукруге? Запишите этот факт.

Повторяйте этот процесс много раз - и вы обнаружите, что отношение количества точек внутри полукруг к общему числу выбрано, назовите это отношение х.

Так как площадь квадрата равна r раз r, вы можете вывести, что площадь полукруга равна x раз r раз r (т.е. x times r squared). Следовательно, x раз 4 даст вам пи.

Это не быстрый способ использования. Но это хороший пример метода Монте-Карло. И если вы посмотрите вокруг, вы обнаружите, что многие проблемы, не связанные с вашими вычислительными навыками, могут быть решены с помощью таких методов.

В интересах полноты, шаблонная версия 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 итераций существует около 80 тыс. Вызовов value() (при отсутствии запоминания).

На самом деле есть целая книга, посвященная (среди прочего) быстрым методам вычисления \pi: 'Pi и AGM' Джонатаном и Питером Борвейном (доступно на Amazon).

Я немного изучил AGM и связанные с ним алгоритмы: он довольно интересный (хотя иногда нетривиальный).

Обратите внимание, что для реализации большинства современных алгоритмов для вычисления \pi вам понадобится библиотека арифметических вычислений с несколькими точными значениями (GMP - неплохой выбор, хотя с тех пор, как я в последний раз использовал его, это было довольно давно).

Временная сложность лучших алгоритмов находится в O (M (n) log (n)), где M (n) - временная сложность для умножения двух n-битных целых чисел (M (n) = O (n). log (n) log (log (n))) с использованием алгоритмов на основе FFT, которые обычно необходимы при вычислении цифр \pi, и такой алгоритм реализован в GMP).

Обратите внимание, что хотя математика, лежащая в основе алгоритмов, может и не быть тривиальной, сами алгоритмы обычно представляют собой несколько строк псевдокода, и их реализация обычно очень проста (если вы решили не писать свою собственную арифметику множественности точности :-)).

Ниже приводится точный ответ, как сделать это максимально быстрым способом - с наименьшими вычислительными затратами. Даже если вам не нравится ответ, вы должны признать, что это действительно самый быстрый способ получить значение PI.

Самый быстрый способ получить значение Pi:

1) выберите свой любимый язык программирования 2) загрузите его библиотеку математики 3) и обнаружите, что Pi там уже определен - готов к использованию!

Если у вас нет математической библиотеки под рукой..

ВТОРОЙ БЫСТРЫЙ способ (более универсальное решение) - это:

найдите Пи в Интернете, например, здесь:

http://www.eveandersson.com/pi/digits/1000000 (1 миллион цифр.. какова ваша точность с плавающей запятой?)

или здесь:

http://3.141592653589793238462643383279502884197169399375105820974944592.com/

или здесь:

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

Очень быстро найти нужные цифры для любой точной арифметики, которую вы хотели бы использовать, и определив константу, вы можете быть уверены, что не тратите драгоценное время процессора.

Это не только отчасти юмористический ответ, но на самом деле, если бы кто-то пошел дальше и вычислил значение Pi в реальном приложении... это было бы довольно большой тратой процессорного времени, не так ли? По крайней мере, я не вижу реального приложения для попытки пересчитать это.

Уважаемый Модератор: обратите внимание, что ОП спросил: "Самый быстрый способ получить значение ПИ"

BBP formula позволяет вычислить n-й разряд - в базе 2 (или 16) - без необходимости даже беспокоиться о предыдущем n -1 цифр:)

Вместо определения pi как константы, я всегда использую acos(-1).

Просто наткнулся на этот, который должен быть здесь для полноты:

рассчитать PI в Piet

Он обладает довольно приятным свойством, что точность может быть улучшена, делая программу более крупной.

Здесь некоторое представление о самом языке

Если эта статья истинна, тогда алгоритм, который Bellard может быть одним из самых быстрых. Он создал pi до 2.7 TRILLION цифр с помощью ПК DESKTOP!

... и он опубликовал здесь

Хорошая работа Беллард, ты пионер!

http://www.theregister.co.uk/2010/01/06/very_long_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 "Costant 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

Вот фрагмент исходного кода для метода 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;
}

РЕДАКТИРОВАТЬ: У меня были некоторые проблемы с вырезать и вставить и идентификацию, в любом случае вы можете найти источник здесь.

Если вы быстрее всего набираете код, то golfscript:

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

Используйте формулу, подобную 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.

Реализовано в схеме, например:

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

Если вы хотите использовать приближение, 355 / 113 хорош для 6 десятичных цифр и имеет дополнительное преимущество использования целочисленных выражений. Это не так важно в наши дни, так как "математический сопроцессор с плавающей запятой" перестает иметь какое-либо значение, но это было очень важно один раз.

С удвоениями:

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

Это будет точностью до 14 знаков после запятой, чтобы заполнить двойной (погрешность, вероятно, связана с тем, что остальные десятичные знаки в дуговых касательных усечены).

Также Seth, это 3.14159265358979323846 3, а не 64.

Пи ровно 3! [Профессор Фринк (Симпсоны)]

Шутка, но здесь один в С# (требуется .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;
        }
    }
}

Рассчитать PI во время компиляции с D.

(Скопировано из 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 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))) );

В этой версии (в Delphi) нет ничего особенного, но она, по крайней мере, быстрее, чем версия, которую Ник Ходж опубликовал в своем блоге :). На моей машине, она занимает около 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.

Если вы хотите вычислить приближение значения π (по какой-то причине), вы должны попробовать алгоритм бинарного извлечения. Улучшение

Если вы хотите получить приближение значения π для выполнения вычислений, тогда:

PI = 3.141592654

Конечно, это только приближение, а не совсем точное. Это немного больше, чем 0,00000000004102. (четыре десяти триллиона, около ⁴/ _{10 000 000 000}).

Если вы хотите сделать математику с π, тогда возьмите себе карандаш и бумагу или пакет компьютерной алгебры и используйте π точное значение, π.

Если вам действительно нужна формула, это весело:

π = -i ln (-1)

В старые времена, с небольшими размерами слов и медленными или несуществующими операциями с плавающей запятой, мы делали такие вещи следующим образом:

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

Для приложений, которые не требуют большой точности (например, для видеоигр), это очень быстро и достаточно точно.

Метод Brent, вышедший выше Крисом, очень хорош; Брент обычно является гигантом в области арифметики произвольной точности.

Если все, что вы хотите, это N-й разряд, знаменитый формула BBP полезен в hex

Один из способов Формула Пи > :

#include <stdio.h>

int main() {
    int r[2800 + 1];
    int i, k;
    int b, d;
    int c = 0;

    for (i = 0; i < 2800; i++) {
        r[i] = 2000;
    }

    for (k = 2800; k > 0; k -= 14) {
        d = 0;

        i = k;
        for (;;) {
            d += r[i] * 10000;
            b = 2 * i - 1;

            r[i] = d % b;
            d /= b;
            i--;
            if (i == 0) break;
            d *= i;
        }
        printf("%.4d", c + d / 10000);
        c = d % 10000;
    }

    return 0;
}

Вычисление π из области круга: -)

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

Лучший подход

Чтобы получить выходные данные стандартных констант, таких как pi или стандартные понятия, мы должны сначала использовать встроенные методы, доступные на языке, который вы используете. Он вернет значение самым быстрым и лучшим способом. Я использую Python, чтобы получить самый быстрый способ получить значение пи

• Пи переменная математической библиотеки. В библиотеке Math переменная 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

• Используйте метод аркос, потому что математика

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

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

Поэтому лучше всего использовать встроенный метод, предоставляемый языком, потому что он самый быстрый и лучший для получения результата. В питоне используйте math.pi