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

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

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

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

Эта версия (в Дельфи) является ничем специальным, но это, по крайней мере, быстрее, чем версия 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.

Пи равняется точно 3! [Профессор Frink (Симпсоны)]

Шутка, но вот один в C# (Требуемая Платформа.NET).

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

Это - "классический" метод, очень легкий реализовать. Эта реализация, в 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)

можно найти больше информации здесь .

Так или иначе самый быстрый способ получить точное as-much-as-you-want значение пи в Python:

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

вот часть источника для gmpy метода пи, я не думаю, что код так же полезен как комментарий в этом случае:

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

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

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

вычисляют PI в Piet

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

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

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

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

Существует на самом деле целая выделенная книга (среди других вещей) к быстра методы для вычисления \pi: 'Пи и AGM', Jonathan и Peter Borwein ( доступный на Amazon).

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

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

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

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

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

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

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

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

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

, Так как областью квадрата являются r времена r, можно вывести, что областью полу круга являются x времена r времена r (то есть, x времена r в квадрате). Следовательно x времена 4 даст Вам пи.

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

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

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

Если самым быстрым Вы означаете самый быстрый вводить в коде, вот решение golfscript :

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

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

Если все, что вам нужно, это N-я цифра, знаменитая Формула BBP полезен в шестнадцатеричном формате

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

Если вы хотите вычислить приближение значения π (по какой-то причине), вам следует попробовать алгоритм двоичного извлечения. улучшение BBP Беллардом дает ИП в O (N ^ 2).

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

PI = 3.141592654

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

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

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

π = - i ln (-1)

В интересах полноты, версия шаблона на 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() (при отсутствии memoisation).

Если эта статья верна, то алгоритм , созданный Беллардом , может быть одним из самых быстрых из доступных. Он создал число Пи в 2,7 ТРИЛЛИОНА цифр с помощью НАСТОЛЬНОГО ПК!

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

Хорошая работа, Беллард, ты первопроходец!

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

С двойными:

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

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

Также Сет, это 3,14159265358979323846 3 , а не 64.