BBP-алгоритм в Python — работа с арифметикой произвольной точности

Пишу курсовую по вычислению числа пи. Пока я закончил теоретический сайт, сейчас я пытаюсь реализовать алгоритм BBP на Python.

Вы можете найти алгоритм BBP здесь: http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula

И это моя реализация на Python:

from sympy.mpmath import *
pi = mpf(0)
mp.dps = 30;
for n in range(0,500):
    pi = pi + mpf((1/16**n)*(4/(8*n+1)-  2/(8*n+4)-  1/(8*n+5)-  1/(8*n+6)) )


print(pi)

Моя проблема в том, что независимо от того, как высоко я устанавливаю k или как высоко я устанавливаю десятичные разряды для числа пи, я не могу получить более высокую точность, чем 16 цифр.

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

Как я могу улучшить свой код, чтобы получить больше цифр?


python algorithm pi arbitrary-precision
person Tim    schedule 01.03.2015    source источник
comment
Мне кажется, что 16**500 просто не может быть представлено float или double, которые вы слишком поздно конвертируете в mpf.   -  person Willem Van Onsem    schedule 01.03.2015
comment
@CommuSoft Что вы имеете в виду, говоря, что конвертируете в MPF слишком поздно? Что я должен делать ?   -  person Tim    schedule 01.03.2015
comment
16**500 не может быть представлено как float, однако, выполнив mpf(16)**500, оно может быть представлено.   -  person Willem Van Onsem    schedule 01.03.2015

Ответы (1)


arrow_upward
1
arrow_downward

По умолчанию python будет использовать стандартные числа с плавающей запятой, определенные IEEE-754. Это имеет точность около 12 цифр и может представлять числа от 2-1022, теперь вы можете решить эту проблему, вызвав оператор mpf ранее в процессе, более читаемая и более точная версия выглядит следующим образом:

from sympy.mpmath import *
pi = mpf(0)
mp.dps = 30;
for n in range(0,500):
    u = 4.0/(8*n+1)-2.0/(8*n+4)-1.0/(8*n+5)-1.0/(8*n+6)
    u = mpf(u)
    d = mpf(16.0/1.0)**n
    pi += u/d

print(pi)

Однако у этого все еще есть проблемы, когда вы вычисляете часть u. Чтобы сделать это точно, вы можете использовать:

from sympy.mpmath import *
pi = mpf(0)
mp.dps = 50;
for n in range(0,50000):
    u = mpf(4)/mpf(8*n+1)-mpf(2)/mpf(8*n+4)-mpf(1)/mpf(8*n+5)-mpf(1)/mpf(8*n+6)
    d = mpf(16.0)**n
    pi += u/d

print(pi)

Что верно для первых 50 цифр:

3.1415926535 8979323846 2643383279 502884197 1693993751

(добавлены пробелы)

person Willem Van Onsem    schedule 01.03.2015
comment
Прежде всего Спасибо! Возможна ли дальнейшая оптимизация или улучшение вычислений возможно только за счет увеличения используемых десятичных разрядов или n ? - person Tim; 01.03.2015
comment
Да, возможно, вы также можете ограничить количество цифр конструкторов mpf в части u или вывести конструкторы констант из цикла. - person Willem Van Onsem; 01.03.2015