Как вычислить число пи c

Как вычислить число пи c

Самый простой и легкий в реализации метод.

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

Будем случайным образом выбирать точки в этом квадрате и считать количество точек, попавших в четверть круга. Благодаря теории вероятности мы знаем, что отношение попаданий в четверть круга к попаданиям ‘в молоко’ равно отношению площадей — пи/4. Вот, собственно, и весь алгоритм. Чем больше взятых наугад точек мы проверим, тем точнее будет отношение площадей.

Вот простенькая программа на Паскале, считающая пи этим способом. Четыре первых знака требуют на моем PentiumII-300 около 5 минут.

uses Crt; const n=10000000; r=46340; r2=46340*46340; var i,pass : LongInt; x,y : real; begin WriteLn('Поехали!'); Randomize; pass:=0; for i:=1 to n do begin x:=Random(r+1); y:=Random(r+1); if ( x*x+y*y < r2 ) then INC(pass); end; TextColor(GREEN); WriteLn('Число ПИ получилось равным: ',(pass/i*4):0:5); ReadLn; end.

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

Слава Богу, этого и не требуется. Поскольку мы хотим найти не точное значение PI, а лишь его приближение с пятью верными десятичными знаками, нам достаточно просуммировать такое количество первых членов ряда, чтобы сумма всех оставшихся членов не превышала 10 -5 .

Остался, правда, открытым вопрос о том, сколько же все-таки членов ряда нужно просуммировать, чтобы получить результат с требуемой точностью?

Ответ на этот вопрос в 'общем виде' выходит далеко за рамки настоящего обсуждения. Это отдельная тема в курсах математического анализа и численных методов.

К счастью, данный конкретный ряд позволяет найти очень простое правило, позволяющее определить момент, когда следует прекратить суммирование. Дело в том, что ряд Грегори является знакопеременным и сходится равномерно (хотя и медленнее, чем хотелось бы). Это означает, что для любого нечетного n , сумма первых n членов ряда всегда дает верхнюю оценку для PI, а сумма n +1 первых членов ряда - нижнюю.

Значит, как только разница между верхней и нижней оценками окажется меньше, чем требуемая точность, можно смело прекращать вычисления и быть уверенным, что как та, так и другая оценки отличаются от истинного значения PI не более, чем на 10 -5 . В качестве окончательного результата разумно взять среднее значение между полученными верхней и нижней оценками. Таким образом, можно предложить алгоритм, приведенный ниже.

Положить n=0, S1 = 0 и S2 = 0; ( начальные установки ) Увеличить n на 1; ( n становится нечетным ) Вычислить S1 = S2 + 4/(2n-1); (S1 - есть верхняя оценка ) Увеличить n на 1; ( n становится четным ) Вычислить S2 = S1 - 4/(2n-1); (S2 - есть нижняя оценка) Если S1 - S2 >= 10^(-5) перейти к шагу 2; ( нужная точность еще не достигнута ) Напечатать результат (S1 + S2) / 2

При реализации этого алгоритма на машине следует помнить, что ряд Грегори сходится достаточно медленно, и поэтому n может принимать довольно большие значения.

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

p	= 12arctan	1	+ 8arctan	1	- 5arctan	1
4		18		57		239

Доказательство этой формулы несложное, поэтому мы его опустим.

Исходник программы, включающий в себя 'длинную арифметику'

Программа вычисляет NbDigits первых цифр числа Пи. Функция вычисления arctan названа arccot, так как arctan(1/p) = arccot(p), но расчет происходит по формуле Тейлора именно для арктангенса, а именно arctan(x) = x - x 3 /3 + x 5 /5 - . x=1/p, значит arccot(x) = 1/p - 1 / p 3 / 3 + . Вычисления происходят рекурсивно: предыдущий элемент суммы делится и дает следующий.

/* ** Pascal Sebah : September 1999 ** ** Subject: ** ** A very easy program to compute Pi with many digits. ** No optimisations, no tricks, just a basic program to learn how ** to compute in multiprecision. ** ** Formulae: ** ** Pi/4 = arctan(1/2)+arctan(1/3) (Hutton 1) ** Pi/4 = 2*arctan(1/3)+arctan(1/7) (Hutton 2) ** Pi/4 = 4*arctan(1/5)-arctan(1/239) (Machin) ** Pi/4 = 12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss) ** ** with arctan(x) = x - x^3/3 + x^5/5 - . ** ** The Lehmer's measure is the sum of the inverse of the decimal ** logarithm of the pk in the arctan(1/pk). The more the measure ** is small, the more the formula is efficient. ** For example, with Machin's formula: ** ** E = 1/log10(5)+1/log10(239) = 1.852 ** ** Data: ** ** A big real (or multiprecision real) is defined in base B as: ** X = x(0) + x(1)/B^1 + . + x(n-1)/B^(n-1) ** where 0 Work with double instead of long and the base B can ** be choosen as 10^8 ** => During the iterations the numbers you add are smaller ** and smaller, take this in account in the +, *, / ** => In the division of y=x/d, you may precompute 1/d and ** avoid multiplications in the loop (only with doubles) ** => MaxDiv may be increased to more than 3000 with doubles ** => . */ #include #include #include #include long B=10000; /* Working base */ long LB=4; /* Log10(base) */ long MaxDiv=450; /* about sqrt(2^31/B) */ /* ** Set the big real x to the small integer Integer */ void SetToInteger (long n, long *x, long Integer) < long i; for (i=1; i/* ** Is the big real x equal to zero ? */ long IsZero (long n, long *x) < long i; for (i=0; i/* ** Addition of big reals : x += y ** Like school addition with carry management */ void Add (long n, long *x, long *y) < long carry=0, i; for (i=n-1; i>=0; i--) < x[i] += y[i]+carry; if (x[i]> > /* ** Substraction of big reals : x -= y ** Like school substraction with carry management ** x must be greater than y */ void Sub (long n, long *x, long *y) < long i; for (i=n-1; i>=0; i--) < x[i] -= y[i]; if (x[i]<0) < if (i) < x[i] += B; x[i-1]--; >> > > /* ** Multiplication of the big real x by the integer q ** x = x*q. ** Like school multiplication with carry management */ void Mul (long n, long *x, long q) < long carry=0, xi, i; for (i=n-1; i>=0; i--) < xi = x[i]*q; xi += carry; if (xi>=B) < carry = xi/B; xi -= (carry*B); >else carry = 0; x[i] = xi; > > /* ** Division of the big real x by the integer d ** The result is y=x/d. ** Like school division with carry management ** d is limited to MaxDiv*MaxDiv. */ void Div (long n, long *x, long d, long *y) < long carry=0, xi, q, i; for (i=0; i> /* ** Find the arc cotangent of the integer p (that is arctan (1/p)) ** Result in the big real x (size n) ** buf1 and buf2 are two buffers of size n */ void arccot (long p, long n, long *x, long *buf1, long *buf2) < long p2=p*p, k=3, sign=0; long *uk=buf1, *vk=buf2; SetToInteger (n, x, 0); SetToInteger (n, uk, 1); /* uk = 1/p */ Div (n, uk, p, uk); Add (n, x, uk); /* x = uk */ while (!IsZero(n, uk)) < if (p/* One step for small p */ else < Div (n, uk, p, uk); /* Two steps for large p (see division) */ Div (n, uk, p, uk); > /* uk = u(k-1)/(p^2) */ Div (n, uk, k, vk); /* vk = uk/k */ if (sign) Add (n, x, vk); /* x = x+vk */ else Sub (n, x, vk); /* x = x-vk */ k+=2; sign = 1-sign; > > /* ** Print the big real x */ void Print (long n, long *x) < long i; printf ("%d.", x[0]); for (i=1; iprintf ("\n"); > /* ** Computation of the constant Pi with arctan relations */ void main () < clock_t endclock, startclock; long NbDigits=10000, NbArctan; long p[10], m[10]; long size=1+NbDigits/LB, i; long *Pi = (long *)malloc(size*sizeof(long)); long *arctan = (long *)malloc(size*sizeof(long)); long *buffer1 = (long *)malloc(size*sizeof(long)); long *buffer2 = (long *)malloc(size*sizeof(long)); startclock = clock(); /* ** Formula used: ** ** Pi/4 = 12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss) */ NbArctan = 3; m[0] = 12; m[1] = 8; m[2] = -5; p[0] = 18; p[1] = 57; p[2] = 239; SetToInteger (size, Pi, 0); /* ** Computation of Pi/4 = Sum(i) [m[i]*arctan(1/p[i])] */ for (i=0; i0) Add (size, Pi, arctan); else Sub (size, Pi, arctan); > Mul (size, Pi, 4); endclock = clock (); Print (size, Pi); /* Print out of Pi */ printf ("Computation time is : %9.2f seconds\n", (float)(endclock-startclock)/(float)CLOCKS_PER_SEC ); free (Pi); free (arctan); free (buffer1); free (buffer2); >

Конечно, это не самые эффективные способы вычисления числа пи. Существует еще громадное количество формул. Например, формула Чудновского (Chudnovsky), разновидности которой используются в Maple. Однако в обычной практике программирования формулы Гаусса вполне хватает, поэтому эти методы не будут описываться в статье. Вряд ли кто-то хочет вычислять миллиарды знаков пи, для которых сложная формула дает большое увеличение скорости.

Глава 6. Вычисление числа π

Одиннадцать первых цифр числа π = 3,1415926535… легко запомнить с помощью такой мнемоники:

«Кто и шутя и скоро пожелаетъ пи узнать число, ужъ знаетъ».

Её придумали до реформы русской орфографии 1918 года, поэтому и буквы «ѥръ» в конце слов после согласных.

Для приближённого вычисления числа π применяются разные методы.

Один из способов (на практике не очень удобный) связан с вычислением бесконечной суммы (ряда) Лейбница: 1 − 1 3 + 1 5 − 1 7 + … = ∑ i = 0 ∞ − 1 i 2 ⁢ i + 1 = π 4 .

Можно получить π из суммы другого ряда — ряда Эйлера: 1 + 1 2 2 + 1 3 2 + … = ∑ i = 1 ∞ 1 i 2 = π 2 6 .

Естественно, на компьютере невозможно осуществить вычисление бесконечной суммы. Однако в математическом анализе доказывается, что оба ряда сходятся. Это значит, что последовательности частных сумм стремятся к некоторым числам (в наших примерах к π 4 и π 2 6 ), то есть сумма достаточно большого количества слагаемых ряда даст хорошее приближение для суммы. Задавшись требуемой точностью, можно даже указать необходимые количества слагаемых, чтобы обеспечить вычисление с этой точностью (конечно, без учёта ошибок округления при выполнении арифметических действий — такие ошибки очень трудно контролировать).

Рассмотрим геометрическую прогрессию со знаменателем − 1 : 1 − 1 + 1 − 1 + … Найдём её сумму S : S = 1 − 1 + 1 − 1 + … = 1 − 1 − 1 + 1 − … . Заметим, что выражение в скобках ничем не отличается от выражения для S : S = 1 − S . Отсюда находим, что S = 1 2 .

Всё это очень странно…

Сходится не всякий ряд. Единственный ряд, изучаемый в школе — геометрическая прогрессия: 1 + q + q 2 + q 3 + q 4 + … = ∑ i = 0 ∞ q i . Его сходимость зависит от q : ряд сходится лишь при q < 1 . В этом случае сумма, как известно, равна 1 1 − q .

Для сходимости необходимо, чтобы последовательность слагаемых стремилась к нулю. Но это условие не является достаточным.

Наряду с бесконечными суммами (рядами) рассматривают также бесконечные произведения и другие бесконечные выражения.

Как вычислить число пи c

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

13 лет назад , # |

В яве длинка с базой равной степени двойки. Если делать с базой 10 k то возможно на вывод потребуется гораздо меньше времени.

13 лет назад , # ^ |

Только мне кажется, что вывести 20000 знаков, преобразовав их в десятичную запись - дело
13 лет назад , # ^ |
Да, только тебе)
Обычный алгоритм смены системы счисления (без быстрого умножения) работает за O(n^2)
13 лет назад , # ^ |

← Rev. 2 → 0

На самом деле нормально написанное преобразование 20000 знаков в десятичную запись за O(n 2 ) работает быстрее, чем за 0.1 секунды. У меня на ноутбуке число из 66666 двоичных цифр переводится в число из 20069 десятичных цифр примерно за 40 миллисекунд. Но это на C++. Попробовал то же самое на Java. Работает 12-13 миллисекунд. Перевод числа из 666666 двоичных цифр в число из 200687 десятичных на C++ проработал 4 секунды, а на Java 1.2 секунды (действительно O(n 2 ) ^^). Хмм. Некоторые программы, оказывается, можно ускорить, переписав их с C++ на Java. 🙂 Хотя, дело скорее всего в 64-битном компиляторе для Java. Если скомпилировать код на C++ 64-битным компилятором под линуксом, то он работает 1.3 секунды.

13 лет назад , # ^ |

>Хмм. Некоторые программы, оказывается, можно ускорить, переписав их с C++ на Java. 🙂

У явы очень мощный JIT оптимизатор, но он помогает только если программа работает достаточно долго чтобы он успел сработать)

13 лет назад , # |

Попробуйте использовать ту хитрую формулу)) Она позволяет находить цифры числа π независимо друг от друга и весьма точно. Идея такая - домножим эту формулу на 16 t и потом целую часть отбросим, а потом в даблах вычислим дробную часть и округлим.

Подробностей не помню - давно решал подобную задачу. В той задаче нужно было вывести i -ю цифру двоичной записи π . Можно ли хитрую формулу эффективно распространить на много цифр - не знаю.

13 лет назад , # |
Топ решений - на Haskell, Python и LISP. Как они это делают??
13 лет назад , # |

← Rev. 3 → +3

I think this algorithm can be useful.
13 лет назад , # ^ |

This algorithm is useful in some way, but, despite very fast increasing number of correct digits, it's slower than the best I tried, 3000 digits against 7000.

13 лет назад , # |

Вообще, при правильном использовании третья("забавная") формула позволяет считать миллионы знаков:

13 лет назад , # |

← Rev. 2 → 0

If you want to try a funny method for calculating Pi digits, I suggest you try some fixed point iterations:
x[0] = 3.14 (or anything close to Pi) and then:
x[i+1] = x[i] + sin(x[i]) -- cubic convergence (i.e. x[i+1] has ~3 times more exact digits than x[i]),
or better:
x[i+1] = x[i] + sin(x[i]) + 1/6*sin^3(x[i]) -- order 5 convergence
or perhaps even:
x[i+1] = x[i] + sin(x[i]) + 1/6*sin^3(x[i]) + 3/40*sin^5(x[i]) -- order 7 convergence.

An advantage of these methods is that they are 'fixed point iterations'. As a consequence, you don't actually need to be precisely accurate in your computations, because the iteration will converge to Pi for any decent approximation. So you should gradually (according to convergence order) increase arithmetic accuracy as the iteration proceeds.
Also notice that only one sin() has to be computed at each iteration, the rest is just a polynomial combination.

So the whole difficulty is now in calculating sin(x). Note that using some standard series around x=0 would be very inefficient, because x will be closer to Pi each time. A funny recursive method to calculate sin(x) is

sin(x, depth) = 3t-4t^3, where t=sin(x/3, depth-1)
sin(x, 0) = x

More depth - more accuracy. Also you could replace base condition by e.g.
sin(x,0) = x - 1/6*x^3,
or include even more terms, then you'll be able to use smaller depth. You'll need to find an optimal combination.

Март, четырнадцатое. Как вычислить число Пи

Каждый год 14 марта мы публикуем материал, посвященный числу Пи. На этот раз по просьбе N + 1 математик Владимир Потапов рассказал о том, как число Пи можно вычислить математическим путем, с помощью различных — подчас необычных — формул.

Еще в древности люди заметили, что отношение длины окружности к ее диаметру близко к трем, но не точно три, а чуть больше. Причем это отношение не зависит ни от диаметра окружности, ни от места, где она проведена. В те времена это отношение, названное впоследствии числом Пи, не сильно выделялось из множества других чисел, которые можно определить опытным путем. Таких как отношение диагонали квадрата к его стороне или отношение площадей квадрата и равностороннего треугольника с такой же, как у квадрата, стороной.

Отцом числа Пи следует считать Архимеда, которого называют автором удивительных открытий, что отношение Пи не приближенно, а в точности связывает не только диаметр и длину окружности, но и площадь круга и квадрат его радиуса, объем шара и куб его радиуса и даже площадь сферы и квадрат ее радиуса. То есть Архимед доказал известные всем со школы формулы: L = 2πr, S1 = πr2, V = 4/3 x πr3 и S2 = 4πr2.

Архимеду принадлежит также первая не опытная, а теоретическая (методом построения описанных и вписанных в круг многоугольников) оценка числа Пи: 3 + (10/71) < < 3 + (1/7). То есть он нашел первые три десятичные цифры числа Пи = 3,14. что и определило сегодняшнюю дату.

Впоследствии математики поняли, что число Пи связывает объем многомерного шара и степень его радиуса при любой размерности пространства (с рациональным множителем, уже зависящим от размерности: для 2х измерений это 1, для 3х измерений — 4/3). Таким образом, число Пи не изменится даже для исследователей, живущих в пространствах с другим числом измерений.

Однако отношение длины окружности к ее диаметру меняется при искривлении пространства и совпадает с нашей константой только в «плоском» однородном случае, проще говоря в пространстве, для которого справедлива теорема Пифагора. Как утверждает теория относительности, рядом с горизонтом событий черной дыры пространство сильно искривлено. Неужели цивилизация, которой повезло возникнуть в подобном месте, может не подозревать о существовании константы Пи?

Оказывается, число Пи неожиданно возникает просто из натурального ряда чисел. Английский математик Джон Валлис, старший современник Исаака Ньютона, открыл удивительную формулу:

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

Еще более удивительную для непосвященных формулу с участием числа вывел великий математик Леонард Эйлер, бóльшую часть своей долгой научной карьеры проработавший в Петербургской академии наук:

Эта формула была признана «самой красивой теоремой в математике». Здесь e = 2,71828. — константа Эйлера, i = √-1 — мнимая единица и Пи — конечно, наше число Пи. На самом деле формулаЭйлера эквивалентна сразу двум равенствам:

где n! = 1×2 x 3···(n — 1) x n.

Конечно, затруднительно вычислять Пи из этих формул как корень уравнения бесконечной степени. А уравнения конечной степени с целыми коэффициентами, корнем которого было бы число Пи, не существует! Это доказал в конце XIX века немецкий математик Фердинанд фон Линдеман, решив заодно знаменитую античную проблему «квадратуры круга». То есть он показал, что, имея отрезок, равный диаметру круга, невозможно только с помощью циркуля и линейки построить квадрат, площадь которого равна площади круга.

Другая знаменитая формула Эйлера:

уже пригодна для приближенного вычисления числа Пи. И даже более подходит для этой цели, чем формула великого немецкого философа и математика Готфрида Лейбница:

Впоследствии выяснилось, что эту формулу задолго до Лейбница вывел индийский математик и астроном Мадхава. Формула Лейбница на самом деле является частным случаем формулы разложения арктангенса в ряд Тейлора:

при подстановке x = 1. Долгое время наиболее удобным для вычисления приближений числа Пи считалось равенство английского математика Джона Мэчина, который был секретарем Лондонского королевского общества, когда его возглавлял Исаак Ньютон. Вот это равенство Мэчина:

Для вычисления числа Пи по формуле Мэчина нужно сначала вычислить arctg1/5 и arctg1/239 с помощью приведенного выше разложения арктангенса в ряд Тейлора, которое, по-видимому, впервые нашел сам Исаак Ньютон.

Число возникает в математике в самых неожиданных местах. Например, математик Абрахам де Муавр (бежавший в Англию из Франции, где его преследовали как гугенота) обнаружил формулу:

Теперь ее называют формулой Эйлера-Пуассона, или интегралом Гаусса.

Сам Муавр, а также выдающиеся математики Пьер-Симон де Лаплас и Карл Фридрих Гаусс в разной степени общности и строгости доказали, что функция Φ(x) = e-x2/2/√2 (из формулы Эйлера-Пуассона следует, что интеграл от функции Φ по вещественной прямой равен 1) является плотностью нормального, или гауссова, распределения, которое является предельным для средних арифметических последовательности независимых случайных величин.

Это означает, например, если мы будем n серий по m раз подбрасывать монету, вычислять разность между числом выпавших «орлов» и «решек» и записывать результат в таблицу, то при росте n и m построенная по таблице гистограмма будет все больше походить на график функции Φ. Эта теорема служит фундаментом для современной квантовой физики, обеспечивая возможность извлекать из многократных измерений случайных событий строгие закономерности.

Казалось бы, тысячелетняя история исследований позволяет предположить, что мы не упустили ничего важного о числе. Однако в 1997 году, совсем недавно в историческом масштабе, произошла сенсация. Саймон Плафф нашел новое представление для числа в виде ряда:

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